{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Bungee jump case study\n",
    "\n",
    "Copyright 2018 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bungee jumping"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose you want to set the world record for the highest \"bungee dunk\", [as shown in this video](https://www.youtube.com/watch?v=UBf7WC19lpw).  Since the record is 70 m, let's design a jump for 80 m.\n",
    "\n",
    "We'll make the following modeling assumptions:\n",
    "\n",
    "1. Initially the bungee cord hangs from a crane with the attachment point 80 m above a cup of tea.\n",
    "\n",
    "2. Until the cord is fully extended, it applies no force to the jumper.  It turns out this might not be a good assumption; we will revisit it.\n",
    "\n",
    "3. After the cord is fully extended, it obeys [Hooke's Law](https://en.wikipedia.org/wiki/Hooke%27s_law); that is, it applies a force to the jumper proportional to the extension of the cord beyond its resting length.\n",
    "\n",
    "4. The jumper is subject to drag force proportional to the square of their velocity, in the opposite of their direction of motion.\n",
    "\n",
    "Our objective is to choose the length of the cord, `L`, and its spring constant, `k`, so that the jumper falls all the way to the tea cup, but no farther! \n",
    "\n",
    "First I'll create a `Param` object to contain the quantities we'll need:\n",
    "\n",
    "1. Let's assume that the jumper's mass is 75 kg.\n",
    "\n",
    "2. With a terminal velocity of 60 m/s.\n",
    "\n",
    "3. The length of the bungee cord is `L = 40 m`.\n",
    "\n",
    "4. The spring constant of the cord is `k = 20 N / m` when the cord is stretched, and 0 when it's compressed.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "newton"
      ],
      "text/latex": [
       "$\\mathrm{newton}$"
      ],
      "text/plain": [
       "<Unit('newton')>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = UNITS.meter\n",
    "s = UNITS.second\n",
    "kg = UNITS.kilogram\n",
    "N = UNITS.newton"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>y_attach</th>\n",
       "      <td>80 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>v_init</th>\n",
       "      <td>0.0 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>g</th>\n",
       "      <td>9.8 meter / second ** 2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mass</th>\n",
       "      <td>75 kilogram</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>area</th>\n",
       "      <td>1 meter ** 2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rho</th>\n",
       "      <td>1.2 kilogram / meter ** 3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>v_term</th>\n",
       "      <td>60.0 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>L</th>\n",
       "      <td>25 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k</th>\n",
       "      <td>40.0 newton / meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>zero_force</th>\n",
       "      <td>0 newton</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "y_attach                       80 meter\n",
       "v_init               0.0 meter / second\n",
       "g               9.8 meter / second ** 2\n",
       "mass                        75 kilogram\n",
       "area                       1 meter ** 2\n",
       "rho           1.2 kilogram / meter ** 3\n",
       "v_term              60.0 meter / second\n",
       "L                              25 meter\n",
       "k                   40.0 newton / meter\n",
       "zero_force                     0 newton\n",
       "dtype: object"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params = Params(y_attach = 80 * m,\n",
    "                 v_init = 0 * m / s,\n",
    "                 g = 9.8 * m/s**2,\n",
    "                 mass = 75 * kg,\n",
    "                 area = 1 * m**2,\n",
    "                 rho = 1.2 * kg/m**3,\n",
    "                 v_term = 60 * m / s,\n",
    "                 L = 25 * m,\n",
    "                 k = 40 * N / m,\n",
    "                 zero_force = 0 * N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now here's a version of `make_system` that takes a `Params` object as a parameter.\n",
    "\n",
    "`make_system` uses the given value of `v_term` to compute the drag coefficient `C_d`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_system(params):\n",
    "    \"\"\"Makes a System object for the given params.\n",
    "    \n",
    "    params: Params object\n",
    "    \n",
    "    returns: System object\n",
    "    \"\"\"\n",
    "    area, mass = params.area, params.mass\n",
    "    g, rho = params.g, params.rho\n",
    "    v_init, v_term = params.v_init, params.v_term\n",
    "    y_attach = params.y_attach\n",
    "    \n",
    "    C_d = 2 * mass * g / (rho * area * v_term**2)\n",
    "    init = State(y=y_attach, v=v_init)\n",
    "    t_end = 20 * s\n",
    "\n",
    "    return System(params, C_d=C_d, \n",
    "                  init=init, t_end=t_end)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's make a `System`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>y_attach</th>\n",
       "      <td>80 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>v_init</th>\n",
       "      <td>0.0 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>g</th>\n",
       "      <td>9.8 meter / second ** 2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mass</th>\n",
       "      <td>75 kilogram</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>area</th>\n",
       "      <td>1 meter ** 2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rho</th>\n",
       "      <td>1.2 kilogram / meter ** 3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>v_term</th>\n",
       "      <td>60.0 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>L</th>\n",
       "      <td>25 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k</th>\n",
       "      <td>40.0 newton / meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>zero_force</th>\n",
       "      <td>0 newton</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C_d</th>\n",
       "      <td>0.3402777777777778 dimensionless</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>init</th>\n",
       "      <td>y              80 meter\n",
       "v    0.0 meter / secon...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_end</th>\n",
       "      <td>20 second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "y_attach                                               80 meter\n",
       "v_init                                       0.0 meter / second\n",
       "g                                       9.8 meter / second ** 2\n",
       "mass                                                75 kilogram\n",
       "area                                               1 meter ** 2\n",
       "rho                                   1.2 kilogram / meter ** 3\n",
       "v_term                                      60.0 meter / second\n",
       "L                                                      25 meter\n",
       "k                                           40.0 newton / meter\n",
       "zero_force                                             0 newton\n",
       "C_d                            0.3402777777777778 dimensionless\n",
       "init          y              80 meter\n",
       "v    0.0 meter / secon...\n",
       "t_end                                                 20 second\n",
       "dtype: object"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "system = make_system(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`spring_force` computes the force of the cord on the jumper.\n",
    "\n",
    "If the spring is not extended, it returns `zero_force`, which is either 0 Newtons or 0, depending on whether the `System` object has units.  I did that so the slope function works correctly with and without units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def spring_force(y, system):\n",
    "    \"\"\"Computes the force of the bungee cord on the jumper:\n",
    "    \n",
    "    y: height of the jumper\n",
    "    \n",
    "    Uses these variables from system|\n",
    "    y_attach: height of the attachment point\n",
    "    L: resting length of the cord\n",
    "    k: spring constant of the cord\n",
    "    \n",
    "    returns: force in N\n",
    "    \"\"\"\n",
    "    y_attach, L, k = system.y_attach, system.L, system.k\n",
    "    \n",
    "    distance_fallen = y_attach - y\n",
    "    if distance_fallen <= L:\n",
    "        return system.zero_force\n",
    "    \n",
    "    extension = distance_fallen - L\n",
    "    f_spring = k * extension\n",
    "    return f_spring"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The spring force is 0 until the cord is fully extended.  When it is extended 1 m, the spring force is 40 N. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "0 newton"
      ],
      "text/latex": [
       "$0\\ \\mathrm{newton}$"
      ],
      "text/plain": [
       "0 <Unit('newton')>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spring_force(80*m, system)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "0 newton"
      ],
      "text/latex": [
       "$0\\ \\mathrm{newton}$"
      ],
      "text/plain": [
       "0 <Unit('newton')>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spring_force(55*m, system)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "40.0 newton"
      ],
      "text/latex": [
       "$40.0\\ \\mathrm{newton}$"
      ],
      "text/plain": [
       "40.0 <Unit('newton')>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spring_force(54*m, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`drag_force` computes drag as a function of velocity:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def drag_force(v, system):\n",
    "    \"\"\"Computes drag force in the opposite direction of `v`.\n",
    "    \n",
    "    v: velocity\n",
    "    system: System object\n",
    "\n",
    "    returns: drag force\n",
    "    \"\"\"\n",
    "    rho, C_d, area = system.rho, system.C_d, system.area\n",
    "    \n",
    "    f_drag = -np.sign(v) * rho * v**2 * C_d * area / 2\n",
    "    return f_drag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the drag force at 60 meters per second."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "735.0 kilogram meter/second<sup>2</sup>"
      ],
      "text/latex": [
       "$735.0\\ \\frac{\\mathrm{kilogram} \\cdot \\mathrm{meter}}{\\mathrm{second}^{2}}$"
      ],
      "text/plain": [
       "735.0 <Unit('kilogram * meter / second ** 2')>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v = -60 * m/s\n",
    "f_drag = drag_force(v, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Acceleration due to drag at 60 m/s is approximately g, which confirms that 60 m/s is terminal velocity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "9.8 meter/second<sup>2</sup>"
      ],
      "text/latex": [
       "$9.8\\ \\frac{\\mathrm{meter}}{\\mathrm{second}^{2}}$"
      ],
      "text/plain": [
       "9.8 <Unit('meter / second ** 2')>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a_drag = f_drag / system.mass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now here's the slope function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def slope_func(state, t, system):\n",
    "    \"\"\"Compute derivatives of the state.\n",
    "    \n",
    "    state: position, velocity\n",
    "    t: time\n",
    "    system: System object containing g, rho,\n",
    "            C_d, area, and mass\n",
    "    \n",
    "    returns: derivatives of y and v\n",
    "    \"\"\"\n",
    "    y, v = state\n",
    "    mass, g = system.mass, system.g\n",
    "    \n",
    "    a_drag = drag_force(v, system) / mass\n",
    "    a_spring = spring_force(y, system) / mass\n",
    "    \n",
    "    dvdt = -g + a_drag + a_spring\n",
    "    \n",
    "    return v, dvdt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As always, let's test the slope function with the initial params."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0 <Unit('meter / second')>, -9.8 <Unit('meter / second ** 2')>)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "slope_func(system.init, 0, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then run the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>success</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>message</th>\n",
       "      <td>The solver successfully reached the end of the...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "success                                                 True\n",
       "message    The solver successfully reached the end of the...\n",
       "dtype: object"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results, details = run_ode_solver(system, slope_func)\n",
    "details"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the plot of position as a function of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_position(results):\n",
    "    plot(results.y)\n",
    "    decorate(xlabel='Time (s)',\n",
    "             ylabel='Position (m)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd1hc17X4/e8UqkCABAghIaG6UUO9S1ax7NhxiaM4dhLHJeXGyc/p17k3ubk3Pa/Ty01z2k0cx3EcxzUuclexZHUQSIKNADWaQIDodWbeP84MHmFAgzQzZ8r6PA+PxHBmZulomDV773X2srhcLoQQQohQYzU7ACGEEGIokqCEEEKEJElQQgghQpIkKCGEECHJbnYAV0opFQcsB2oBh8nhCCGEGB0bMBE4oLXu8f5B2CcojOS0y+wghBBCXJH1wJveN0RCgqoFeOSRR8jKyjI7FiGEEKNQV1fHHXfcAe73cm+RkKAcAFlZWUyePNnsWIQQQlyedyzRSJGEEEKIkCQJSgghREiSBCWEECIkSYISQggRkoJaJKGUWgX8L6CABuB7Wus/KKVigV8Ct2IslP1Ea/1AMGMTQggRWoI2glJKWYFngP/VWqcAHwR+qZRaCHwTI2nNwLiu6W6l1F3Bik0IIUToCeYUXxqQCViUUhbABfQDvcDdwHe11s1a61PAj4B7gxVYe1cfJSebaG7rRtqPCCFEaAjaFJ/WulEp9UvgIeBPGNtbfA7j4qyJwHGvw0uBBcGK7WePHmbfsToAkhNjmJI1lk1Lc9i8LIcYuyzTCSGEGYKWoNxTfN3Ah4AngDXAk8AF9yGdXod3AonBim3zshwutPdQda6Nts4+jlU2cqyykb+/orl100yuWTmV2BhbsMIRQghBcIsktgJrtdZfcn+/Qyn1R4zpPYAEr2MTgfZgBbYmP5s1+dm4XC6aWrspKj/P46+d4Oy5Nh58qph/vXmSr35kBTkTkoMVkhBCRL1gzl/lAHGDbuvHqOarwyiS8Mjj4im/oLBYLIxPSWDT0hx+ef8mvnL3cnImJFHd0M79/7uT/cfrgh2SEEJErWAmqJeBeUqpTyilLEqppcC/AY8CDwNfV0qlK6Vygfvdt5nGarWwJj+bH39uA2vzs+ns7uc7/7ePf7xaJoUUQggRBEFLUFrrYxjTfPdirDv9Dfiy1voZ4GvAUeAYcABjjerBYMU2koQ4O/951zI+fH0eAA+/WMI/Xi0zOSohhIh8Qb1QV2v9AvDCELd3A/e5v0KOxWLh9i2KSRlJ/ODhg/x1WynjUxLYsmKK2aEJIUTEkhrqUVi3cBKfuMWofv/l44Uc1vUmRySEEJFLEtQo3bhuOu/bNBOH08X3HtpPZXWL2SEJIUREkgR1Ge5691w2LJ5MV4+DHzx8kO7efrNDEkKIiCMJ6jJYrRY+c/uigRL0h54PekW8EEJEPElQlykuxsYXP7QUm9XCc2+epEDWo4QQwq8kQV2BmZNT+eC1xvXFP3+sgPbOXpMjEkKIyCEJ6grdunkWamoajS3d/PapYrPDEUKIiCEJ6grZbFa++MElxMbY2H64imOVjWaHJIQQEUESlB9kZySxdeNMAP7wTDFOp2yFJIQQV0oSlJ+8b9NMxo2Np7yqhe2Hz5odjhBChD1JUH4SH2fn7hvmAPDQ8yV098i1UUIIcSUkQfnRxiU5zMxJpam1myfeKDc7HCGECGuSoPzIarXw8ZvnA/Dk9nIaW7pMjkgIIcKXJCg/mzd9PKsXTKS3z8GTMooSQojLJgkqADwX727be5oLbT0mRyOEEOFJElQATMtOYcXcLHr7HDyzs8LscIQQIixJggqQ27bMAuD53SdlCyQhhLgMkqACRE0dx8JZ6XT19PPc7pNmhyOEEGFHElQA3b7FWIt6dmcFXXJdlBBCjIokqACaP2M8c3LH0dbZx7a3TpkdjhBChBVJUAFksVi49WpjLeq5NytxyB59QgjhM0lQAbYsbwJZ4xOpb+7iUMk5s8MRQoiwIQkqwKxWC+9eMw0wRlFCCCF8IwkqCLasmEKs3UpBWQPVDe1mhyOEEGFBElQQJCfGsmHJZABe2CMl50II4QtJUEFyw1pjmu+1/WekFYcQQvhAElSQzJicypzccXR097P9cJXZ4QghRMiTBBVE73aPop7ffRKXS0rOhRBiJHazA4gma/Oz+f3TxZyqbaWiqoWZOalmhyRCSE+fg5KTjRSWNXCuqZPE+BgS4+0kJ8aybM4Epk9KMTtEIYJKElQQxditbFwymWd3VfLawTOSoAQAJ2taePjFEgrLGujrdw55zMMvljA9O4WrV+SweWkOSYmxQY5SiOCTBBVkm5fl8OyuSnYcruajN80nxi6zrNGqvauPR7aV8MLuk3g2GZk+KYVFszLIzR5LT6+Dzu4+6ho72VVYTWVNC5VPt/D4qyf4zO2LWDE3y9x/gBABJgkqyKZPSiF34lhO1bZysKSO1QuyzQ5JmOBwaT0/efQQLe29WC1w49pp3HbNbNKS44c8/uPvmc++Y3X8a1clJaea+PYf93H9mlw+etM84mPl11hEJvn4HmQWi4Wrl+cA8NqBsyZHI8yw/XAV3/rjXlrae5k3fTw/++JG7t2aP2xyAoiNsbF+0SS+d986PnLjPOw2Cy/uOcUXfrqD+ubOIEYvRPBIgjLBhiWTsVotHCw5Jy3ho8y/dlXy40cO4XC62LpxJg/8v7VMy/a9+MFqtbB100x+8vkN5ExIpqq+nf/+zR4aW7oCGLUQ5pAEZYK05HiW5mXicLrYUSDXREWLx17V/O7pYgA+cuNcPnLTPCwWy2U91rTsFH7wmfXMnJxCbWMHX/3Nbppbu/0ZrhCmkwRlkquXTQHgdZnmiwq7Cqv564ulWC3wudsXsXXTrCt+zKSEGL517xqmZY+luqGDrz64R0bkIqJIgjLJinkTSEqIobKmhZM1LWaHIwLoZE0LP3+sAICPvWc+W1ZM9dtjJyfG8u171zA1K5mz59r4wcMHpe+YiBiSoEwSYzcWvcH4dC0iU2tHL9/50356eh1sXpbDTeum+/05UpLi+Pa9a0hLjqO44jyPvaL9/hxCmEESlInWLTJKzN8srJGtjyKQw+nihw8fpL6pk5k5qdx368LLXnO6lLSx8fz7HUuxWODvr2iOnGgIyPMIEUySoEw0b3o6qUlx1DZ2UFEt03yR5l+7Kik80UBqUhz/dfcKYmNsAX2+hbMyuH2LwuWCHz9yiOY2KZoQ4U0SlIlsVgtr8icCsPtIjcnRCH8619TJX7eVAPCZ2xaRkZYQlOf9wLWK+TPG09zWw8/+XiAjcxHWJEGZbJ17HerNI9XyZhIhXC4Xv3niCD29DtYuzGbFvOBtSWSzWrj/jqUkJ8ZwuLSeN+WDjwhjkqBMNnfaeFKT46hr7KSiSqb5IsGuwmoOldYzJt7OvbcsCPrzj09J4O4b5gLwh2eO0tndF/QYhPAHSVAms1ktrM13F0sckWq+cNfW2cvvnz4KwEdumkfa2OG3Lwqka1ZMZfaUVJpau3nslTJTYhDiSkmCCgHrFhoJatcRqeYLd49sK+VCew/zpo/nGj9e7zRaVquFT27Nx2KBZ3ZWcKau1bRYhLhcQU1QSqmJSqmnlVItSqlzSqlvu2+PVUr9TinVpJRqUEp9JZhxmW3OtPGMGxtHfVMn5VUXzA5HXKa6xg5e2nsKiwU+tTUfqzUwJeW+mpWTxrtW5eJwuvjtU8Xy4UeEnWCPoJ4BaoEJwCrgbqXUh4BvAgqYASx3335XkGMzjc1qYc2Ct6+JEuHp0Zc1/Q4XG5dMZurEsWaHA8Cd188hOTGWovLz7D1aZ3Y4IkR09fRTcrKJQ6XnqKpvo6/fYXZIQwpaIxml1EpgOrBWa90HnFRKbQS6gB8B92itm4FmpdSPgHuBvwQrPrOtWZjNc7tP8tbRWu65cW7ALugUgXG6rpU3Dp3FZrXwoXflmR3OgLFjYvngtYrfPV3Moy+XsnJelukjO2GO03WtPL29gpJTjdSc78B7QG2xQEZaIpuWTObGddNJTY4zL1Avwex0thQoBr6hlLoH6AZ+DfwRmAgc9zq2FAh++ZOJ5uaOIzkxhtrzHVTVt5MzIdnskMQoPLKtFJcLrl09lazxY8wO5yLvWjWVJ944wcmaVt46WjtQlCOiQ01DO397SbOzsGogKdmsFqZMTCY5MZZzTZ00XOiivqmTx14t48nt5WxelsP7r57NhHGJpsYezAQ1DlgP7MAYSeUB2wDPnizeXdc6AXPPTJDZbFaWzZnAG4eq2H+sThJUGCk708xbxbXExtj4wDXK7HDeITbGxvuvns2DTxbx6EulrJ4/UUZRUcDlcvGPV8v428sap9OF3WbhXaty2bJiClOzkomxv72zSb/DiT7dzFPby9l3rI6X9p5mV2E1n//AElYvmGjavyGYCaoHaNVaf8P9/RGl1B+Au93fe19qnwi0BzG2kLBy3kTeOFTFvmN1vG/zlbdjEMHxyLZSAG5aN41xJpWVX8q1K6fwz9dPcLqujd1FNQMbFYvI1NPn4H8fK2BnQTUWC1yzYgofuEaROcyIyG6zMm/6eOZNH8/Zc2089Pxx9h2r4//7835u3TyLD1+Xh80W/KLvYD5jKZColIr1us0ONAN1GEUSHnlcPOUXFRarDOw2K6Wnm2hpl74+4eBkTQuHdT0JcTa/9HgKlBi7jdu3zAbg0ZdLpSVHBGtu7earv97NzoJqEuJs/PdHV/LZ2xcPm5wGy5mQzFc/soKP3DgPqwX++foJvvH7vXT19Ac48ncKZoJ6BWM678fusvIFwMeAR4GHga8rpdKVUrnA/e7bokpifAz5s9JxueDA8XNmhyN88PSOCsC4MHbsmNhLHG2uq5dPITMtgbPn2qXFS4Rqae/hP3/1JvpMMxlpCXz/0+tZMXf0W21ZLBa2bprJtz+5htSkOApPNPC9hw7Q73AGIOrhBS1Baa27gQ0Y60+1GOtPP9BaPwF8DTgKHAMOAE8ADwYrtlCy0r1v2/7jUhIc6ppau9lZUIXFAjet93+fJ3+LsVu5bYsxUfHkGyfkuqgI09vn4Lt/2k/t+Q6mZ6fw489dxbTslCt6zPyZGXzv0+sYOyaWw7qenz9WgDOIo+9grkGhta4Ebhji9m7gPvdXVFsxN4vfPFHEYV1Pb58j4C0axOV7fvdJ+h0uVi+YGHKVe8PZtHQyf32xhJM1rRytaGTBzHSzQxJ+4HK5+PljBZScaiI9JZ6vfXwlacn+WQ+dlJHE1z++iq/+ZjfbD1WRmhTHx26e75fHvhTZ6ijEpKcmMGNyCj29Dmk6F8K6e/t5cc8pAG7ZMMPcYEYhNsbG9WtyAWMLJBEZHnmpdGDN6WsfX8X4FP+2d5k9JY3/umcFdpuFp3dU8MKek359/OFIggpBK91zxvuOyTRfqHrjUBVtnb3MnpLKnNxxZoczKtevycVus7L/eB215zvMDkdcof3H63jslTKsFviPO5df8bTecBarTD57+2IA/vjM0aDs7ygJKgStnG9cd3DgeF1Q53uFb5xOF8+4iyNuuWpm2O36kZYcz4Ylk3C54Lk3K80OR1yB9s5efvV4IQB33zCXZXMmBPT5Ni3NYcvyKfT2O/nxI4fp6w9s0YQkqBA0LXss41PiaWrt4VSt7EIdagrLGqhuaCc9NWGgI3K4uXm9MS35yv4z0i8qjP3u6WKaWnuYkzuO92yYGZTn/Ldb5pM1PpHKmhYecXeNDhRJUCHIYrGwRGUCcKhUys1Dzcv7TwNw3aqpply86A/TJ6WwYEY6XT39vLL/jNnhiMuw/1gdbxyqIjbGxuc/sBhbkHYHSYyP4YsfXIrVAk9uL6e4/HzAnsvn3y6l1ASl1LuVUvcope5USl2jlBofsMii3NI8Y6h+WNebHInw1tLew76jtVgtxnVF4ezmq4zS+H/tqpSp5DDT1tnLL91Te3e9ew7ZGUlBff4508Zx2xaFy2UUaATKiGXmSik78CHg88BCoBdj5wcbxt56KKX2YWz6+netdXCv4opgC2dnYLVaKDnZRGd3H4nxMWaHJIAdh6vod7hYmpdJeqp/K6WCbfncLDLTEjjX1ElReQOLZmeaHZLw0UPPH6e5rYe508Zx0zpzrsG7/ZrZ9DucTA9QUQaMMIJSSm0AioC7MHYcnw0kaq2ztdYTgFhgMfA34NNAqbt9hvCDpIQY1JQ0HE4XR04EbggtfOdyuQamw8zslusvNquFLe5R4Cv7ZJovXJypa+WVfaexWi18+v2LTNv4126zcvcNc1m/OHD7Oo40gvp34HatdfFQP9RauzB2fzgK/FoptRj4FrDd30FGq6V5mZScauKwrjd1R2FhKK+6wKnaVsaOiWXFvNFvHxOKrl4+hUdf0bx1tJa2zl6SE0N7uyYBDz1fgtMF16+eGvFdD4YdQWmtbx4uOQ1zfIHW+ib/hCUAluQZUy6HS8/JtjQhwDN62rh0MjH28CyOGCxzXCILZ2XQ1+9k+6Eqs8MRl3C04jz7j9cRH2vjg9eGXmsXf/N5qyOlVCIwDXhHq0Wt9WF/BiUMMyalMnZMLPXNXdLE0GQ9fQ52HjbewK+NgOk9b9eumEphWQOv7j8TFnsKRiuXy8WfnjsGwNaNM/22lVEo8+ljoFLqw0A9xprUwUFfBwIWXZSzWt8uN5dqPnPtKaqho7uf2VNSmTpxrNnh+NWqBVkkJ8ZQWdNCedUFs8MRw9hdVEPZmQukJsdxy8bgXPNkNl/nKR7AKJSYjtGe3ftL+kcH0NvTfJKgzOSZ/toS5qXlQ4mx29i4NAeAV/adNjkaMRSHw8lfXjAuiv3QtYqEuKDu820aX/+VY4Ffaq3l1Rtki92lv0crztPT5yBOdjcPupb2HgpPNGCzWli7MDI70V6zYgr/2lXJjsNVfPTm+fI6CzF7imqpPd9B1vhErlkZWVPMI/F1BPUwcE8A4xDDSE2OY+bkFHr7nRytkHJzM+wpqsHpdLFodkbINyW8XNOyU5iZk0pHdz97i2vNDkd4cblc/PP1E4Cx9mQP091LLoevI6gfAoeVUncAp4CLLsjVWm/2c1zCy2KVSXlVC4VlDQM7TIjg2VVYA8BVAbzeIxRsXppD+dkL7CyoZsOSyWaHI9wKyhqorGkhNTku7HcvGa3RjKDagecxiiIODfoSAbRwZgYARXLBbtA1tnRxtPI8MXYrK+dF9rVo6xZlY7UY+z+2dvSaHY5we8I9erp5/fSoa2Dq6whqObBSa10UyGDE0PKmjSPGbqWypoWW9h5Skt5R6S8CZPeRGlwu46LpMQmRvd1UWnI8+bMyKCxrYE9RDdetzjU7pKhXdqaZovLzJMTZuX7NNLPDCTpfR1AaSA1kIGJ4cTG2gaZ4RQHcOVi8087CagCuWhQdU14bFhv/zp0F1SZHIoCBtafrV+eSFOEfkIbi6wjqAeDPSqlfAhXARQ1ktNYv+DswcbFFszMoKj/PkRMNrF8U2WshoeJcUyf6dDNxsTaWz42Otb/VCyby6yeOcLTyPOcvdIX9hrjhrLqhnb1Ha7HbrAM7z0cbXxPUo+4/fzTEz1wYu5uLAFo4KwMokXWoIHrTPXpaMTeL+Ci57mRMQgzL505gT1EtuwqreW+UXBAail7YfRKXCzYtncz4lOj8oODTb53WOnrqGkPUjMmpjIm3U9vYwbmmTiaMSzQ7pIjnmd6LthHrhsWT2VNUy46CKklQJunu6ee1A8bejzea1E4jFFyq3caoKKWk3DxAbFYL82ekA3DkRIPJ0US+usYOKqtbSIizsTQvuvokLZszgcR4OxVVLZw912Z2OFFpR0EVHd395E1NY/qkwPVbCnUjjYy+oJR6USl1vVJq2NU5pZRdKXWLUupVjMaGIkAWzTbKzSVBBd7eo3WA0dk42kp7Y2NsrFlg7GAmxRLB53K5eH73SQBuWBt9lXvehp3i01rfopR6L/A9YKpSajtwDDgPWIAMjC67q4EzwLe11v8MeMRRzFiHMir5XC4XFos5jcqiwd6jxm4K0dqHa/3iSbx64Ay7i6q547o8s8OJKqWnmjlZ00pKUixrF0b3VqcjrkFprZ8CnnJ3yn03RjKagLGTRB3GRboPaK13BThOAUzOTGLc2DiaWns4U9cWcbtqh4qW9h5KTjZit1midueO/JnpJCfGcPZcO6frWpmaJa+1YPGMnq5dOZUYe3SN3gfztUhiO9Ip13QWi4X8WRlsP1RF4YkGSVABsv9YHU4XLJqZEfEX5w7HbrOyav5EXtl/hj1HaiRBBcmFth52F1VjtcB1q3LNDsd0Up0XZjzbHsk6VOB41p9WRen0nodneml3UY3JkUSPV/afpt/hYvncLDKlUlcSVLjJn2VU8h2vbMThlDbw/tbV009BmdF7a+W8LJOjMVe+ewR5uq5NqvmCwOl08bK7H5dsM2WQBBVmMtMSyRqfSEd3P5XV0v3U3wp0PX39TtTUNMaNjfyW2iMxNsg1kvSeYhlFBdqxk43UNXaSnhLPYhVdlzYMRxJUGFrgvh6qWPbl87u3PNV786N7es/DM82354j0iAq0V/cbF+ZuXj4Fm1UqdMH3rY5QSmUC+UAMRpn5ANmLL7jyZ6bzyv4zFJWfZ+umWWaHEzH6HU4OHD8HyPqTx+LZGSTG26msaaHmfDvZ6UlmhxSROrv7Btb6rl6eY3I0ocOnBKWU+hjwa4zkNJjsxRdkC2a616FONuJwOLFFUYfNQDpW0UhHVx85E5KYlCFvxAAxdhsr5mWx/VAVu4/U8P6rZ5sdUkR680gNPb0O5s8YLx8CvPj6zvYl4PdAitbaOuhLklOQjU9JIDt9DF09DsqrZB3KX/aXGNV7K+ZGd3HEYGvzpZov0DzTe1uirGPupfiaoHKAn2utpZQnRHhGUdIfyn8Ouqf3lkuCushilUlCnI2KqhbqmzrNDifinD3XRsmpJhLibAMfBoTB1wT1MnB1IAMRo5M/Uwol/Km6oZ2a8x0kJ8aQNzXN7HBCSlyMjSXuHTU8W0AJ//HsWr5u4aSoaeviK1/PxhHgJ0qpm4EyoNf7h1rr//B3YGJknkq+46ea6Hc4scs61BU5cNyY3luiJsia3hBWz5/I7iM1vHW0lpuvmmF2OBHD4XDyxqGzAFyzYqrJ0YQeXxPUBmAfkICxQaw3uVrUBGlj45mcmURVfTsnzlxgzrRxZocU1g4MTO9F5957l7JszgTsNgvHKxtpae8hJSnO7JAiwpHy8zS19pCdPoa8XBm5D+brXnybAh2IGL0FM9Opqm+nqKJBEtQV6Ojq41hlI1arhSVR1vvJV2MSYsifmcFhXc+B43VskU/7frHdPXrauDRHuhMMwee5DKXUBKXUt5VSTyqlnlZKPaCUit5WjyFA1qH8o6CsHofTxZzccSQnxpodTsjyXBv2VnGdyZFEhu6eft4qNtb0Ni6ZbHI0ocmnBKWUWoGx9vRejH5QDcCNQJFSalngwhMj8axDlZxsoq/fYXI04Wtgem+OTO+NZNW8LCwWI6F39fSbHU7Y23usju5eB3lT05iYPsbscEKSryOoHwOPAgu01p/QWv+b1noB8Gfgh4EKTowsJSmOqVnJ9PY7KTsj10NdDqfTxaFSI0Etk/WnEaWNjSdv6jj6+p0cLq03O5yw5z29J4bma4JaBvxUaz24IOIXwHL/hiRGY757FHW0Qqb5LseJs820tPeSOS6RKROSzQ4n5K2a75nmk3LzK9Hc1k1BWQM2q4V1Ud41dyS+JqhaIHeI26cDcvGuiRYMJKhGkyMJT57pvRVzJsgitQ9WLTAuYj5YUkdfv9PkaMLXrsJqnE4XS/MmSEXkCHwtM38Y+J1S6vPAXvdtq4Gfun8mTDJv+njAuB6qr99JjF2u4RmNg+7pvaWy/uST7PQkpmYlc7qujeKK8yyRthCXZfuhKgA2LpXiiJH4+m72XYzdJP4BVAHVGGtSjwNfDUxowhepyXFMyUqmt8/BibPNZocTVprbuqmoaiHWbh3YOkpcmqeaT3aVuDxV9W2cOHuBhDg7K6K8Keal+JSgtNa9Wut/A9IxRk4LgVSt9f1a677RPKFSKlUpdUYpdY/7+1il1O+UUk1KqQal1FdG+W+IegP9oWQdalQKtLHQP39mOnExsuexr1bNMxLU/mN1OKWr86jtLKgGYE3+RHndXcKwCUop9W6lVIzX398NrMVIUjnAJq/bR+NBYJLX998EFDADo+DibqXUXaN8zKg2f4YxzXe0XNahRuNQiZGglsrFuaMyY3IK6SnxNLZ0y276o+RyuQYS1FWLZXrvUkZag3oOyALq3X8fjs/9oJRSdwNjgWKvm+8G7tFaNwPNSqkfAfcCf/HlMQXMn+6+Huq0rEP5yuF0UVDmSVCy/jQaFouFlfMn8vzuk+w7VsfsKbJFj68qq1uobmgnJSmWhTKtfEnDJiittXWov18updQ04OvAGmCb+7ZUYCJw3OvQUmDBlT5fNElNjiNnQjJnz7VRflb25fPFibPNtHX2kTU+kWy5SHLUVs7L4vndJ9l7tJY7r59jdjhhY1ehMXpam58tmxL7wNedJF53J5PBt2copQ75cH8b8Ffgfq219z4pntaR3k1mOoFEX+ISb1vgnuaTdSjfvD29J+Xll2P+jHTGxNs5U9dGzfl2s8MJC06ni52FMr03GsOOoJRSG4G57m83APcqpQZf8zQHY+3oUv4H0FrrJwfd3uH+M8HrtkRAXvGjtGBmOi/sOUVxxXlu2yJtuS/Fs3uEbA57eWLsVpbOmcDOgmr2Ha3jvRtnmh1SyNOnm2lo7iI9JZ45uTLL4YuR1qAagfsBi/vrPsB7wzcXRiL5dx+e5wNAtlJqq/v7ZODXwAqgDqNIotr9szwunvITPvBcD1Ui/aEuqaW9h/KqC9htVvJnyDrA5Vo1f6KRoI5JgvLFzgLj2qf1iydjtcqo3RcjrUEVY+wUgVLqDWCru5Bh1LTWed7fK6UKgZ9prf+slGoHvq6UKsKY8rsf+PnlPE80S0uOJ2dCEmfPtVN+9gJ58gltWAW6HpfLqH6UDqaXb2leJnabhZKT0iPqUhwOJ28eqQHgqkWTLnG08BipzNx7HegGoEcpldOb0zEAACAASURBVDjU1xXG8DXgKHAMOAA8gVGKLkbJsy9fkbTfGNGhUikv94fEeKNHlNP1dkdiMbTiivNcaDcaE86YnGJ2OGFjpI+PbUqpiVrreoypvKGuyLMwijJzD631Iq+/d2NMH943mscQ77RgRjov7jnFUVmHGpbT6eKwlvJyf1k1P4vDup69R6WJ4Ug81z6tXzxJinJGYaQEtRlocv9dOuqGAc8Fu7Iv3/Aqqi/Q2tFLRloCkzOTLn0HMaIV87L49RNFFJQ10N3bT3ysTJkO1tfvZI9793eZ3hudkdagdgz1dzC2JwLygTKtdWvgwhOjMXgdSq6HeifP6GmJypRPsn4wPiUBNSUNfaaZAt3Aavc+feJthWX1dHT1MTUrmSlZY80OJ6z4eh3UTKXUDqXUKvea037312ml1KqARihGRfblG9lhWX/yu5XzjQ1PZfPYoXkuzl0vo6dR83UO6BcYfZ9OAXcCkzFKw38D/CQgkYnL4tmVWxLUO3V09VF6uhmr1UL+zAyzw4kYniaGB47X4XBIjyhvvX0O9h41CkgkQY2erwlqPfAF9y4QtwDPa61PAL8HFo14TxFUA/vyudehxNuOnGjA6XQxJ3ccYxJizA4nYkzOTGJSxhjaOvs4frLp0neIIodK6+nq6Wf6pBSyM2TNc7R8TVDdQIxSagzGrhIvum/PAloCEZi4PJ59+Xp6pT/UYJ71p8VKRk/+ZLFYBkZRMs13sTdleu+K+JqgXsIYLT2BsVfev5RSV7tvezZAsYnLJPvyvZPL5VVerqS83N+8E5TLJT2iALp7+tnnvj5s3cJsk6MJT74mqHuBgxgjqRu01h0YvZu2A58PTGjicnnWoaQ/1Nuq6ttpaO4iJSmW6ZPkQkl/mz0ljdTkOOqbuzhZI4W9AAdKztHT62D2lFSyxsuO+ZfDp4sWtNbtwOcAlFJjlVKpWuvvBTQycdk861ByPdTbBqb3ZmfKPmgBYLVaWDkvi5f2nmbf0Vr5EIBU7/mDz+9cSqlPKaXOAs1Ao1KqVin15cCFJi5XanIcU7KS6e2TdSgPT3m57F4eOG9P88m2R53dfRwqMXbMX5svCepy+Xod1P3A9zDKzdcDVwE/Bf5DKfW5wIUnLtfA9VCyLx89fQ6OutfjFs+WBBUoC2elkxBno7KmhbrGjkvfIYLtP1ZHb7+TObnjyEhLuPQdxJB8HUHdB3xSa/0DrfUerfVurfUPgE8Bnw5ceOJyLZCNYwccq2ikt9/JjMkppCbLjtuBEmO3sWyOXLQLDDQm3LBYRk9XwtcElYGx0/hghzAu2hUhxrMvX+mpJnr7HJc4OrJ5b28kAsuz1dGeouhNUO2dvRToeqwWWCPVe1fE1wR1FHj/ELffDpT6LxzhLylJceROHEtvvxN9OrrXoQ5rd/dcSVABtzQvkxi7ldLTTTS1dpsdjineKq6l3+Fiwcx00pLjzQ4nrPm69fDXgOeVUquBt9y3rQauA7YOey9hqvxZ6ZyqbeVIecNA6Xm0qW/q5Oy5dhLi7NLEMQgS42NYPDuT/cfr2He0luvXTDM7pKDbOVC9J5NLV8qnEZTW+mXgaqAHYy++W4FWYLnW+rnAhSeuxEL3fnNFJ6J3HcozvbdodgZ2m5TbB8PANF9x9E3zXWjroaj8PDarhTX5srP7lfK5eYvWeiewM4CxCD+bN308VguUnWmmq6efhChsby7rT8G3Yl4WVquF4vLztHX2kpwYa3ZIQbOnuAan08WyOROi6t8dKCO2fFdK/U4p1eS+5unXSilpZhJGxiTEMDMnFYfTxfGT0berRL/DSWFZAyAJKpjGjollwYzxOJyuqGsFLxfn+tdIcx7fBG4CfoDRUuMGjL33RBjJj+JpvtJTTXT19JMzIYnMcYlmhxNVVi8wqteiqZqvsaWLY5WNxNitrHL3yBJXZqQEdSvwIa3197TWP8So4nuPUkr6FISR/Jme66EaTI4k+N7evVxGT8HmeYMu0Ea7iWjw5pEaXC5YNmcCifHyNukPIyWoyVxcQn7AfbxsBR1G5kwbh91moaK6hfbOXrPDCSrZvdw841MSUFPT6O13cqj0nNnhBMWOw1UAbFgs1Xv+MlKCsgEDV3hqrV0YVXyy8hdG4mPtqKnjcLngaGX0rEM1t3VTUdVCrN3KPPdFyyK41uYb03y7j9SYHEng1TS0c+LsBRLi7CybKx+I/EXqbqPAwpnRt+1RgTamNOfPTCcuxmZyNNHJk6AOlJyjO8Kn+Tyjp9ULJsrrzY8uVXd8j1KqfdDxH1ZKXfROp7X+td8jE36TPyuDv72sKToRPetQnmmlpbL+ZJrMcYmoqWno080cLD3HuoWRWdnmcrnYUeCe3lsi03v+NFKCOoOxGay3OuAjg25zAZKgQtjsKWnExdo4XddGc2s3aWMje/sVh9NFgZb2GqFg3cJJ6NPNvFlYE7EJqqKqheqGDlKT4gZmK4R/DJugtNa5QYxDBFCM3cq86eM5XFrPkRMNbFyaY3ZIAXXiTDNtnX1kjU9kUkaS2eFEtbX52fzx2aMcKDkXsReLe0ZP6xZlY5PdSvxKzmaUWDzbuB6qoCzyp/kOuhvFLcubgMUi3XPNlJGWwJzccfT2OTh4PPKq+RxOFzsL3K01ZHrP7yRBRYlF7kZ9hWX1uFwuk6MJrIH1pzlSTRUK1rlbTuw6Um1yJP53rPI8Ta3dTBiXiJqSZnY4EUcSVJSYmpVMWnIcTa09nKlrMzucgGlu7abcXV4erTu4h5o17mq+QyXn6OzuMzka/9p+6O3iCBmt+58kqChhsVhYFAXTfIdKjeKI/FkZUu4bItJT3dN8/U4ORNA0X0+fg91FxjVe0jk3MCRBRRHvab5IddAzvSfVeyFl3SL3NF9h5Ezz7S2upbO7n5k5qUzJkn20A0ESVBTxjKCKKxrp64+8NvAOh5NCd3n5Mll/CinrFk7CajHWB1s7ImPLrdcPngXg6mWRXRVrJklQUWTc2HijDXyfg5JTTWaH43elp5vp6O5nUkYSWePHmB2O8DJubDwLZ2XQ73CxOwKKJRpbuigsq8dus3CV7L0XMJKgooxnFFUYgetQA+XlMnoKSZvcI4033IUF4Wz7oSqcLlg+N4uxY2R70kCRBBVlFrvXoTw7LUSStxOUrD+FolXzJxIXa6PkVBN1jR1mh3PZXC4Xr7mn9zbL9F5ASYKKMnOnj8Nus1JR3UJLe4/Z4fhNQ3MXp2pbiY+1MW+67F4eihLi7KxeMBGA7YfDdxRVUdXC2XNtjB0Ty9I8Ga0HkiSoKBMfa2fuNKP9RiR12d3vbi2+WGUSY5fy8lC1aYl7mu/g2bC9YPy1g2cA49qnGLu8hQaSnN0otMS9w/fBCGok50lQK+ZKq+1QtnBWOmnJcdSc7+DE2QtmhzNqff1Odhw2ijxkei/wJEFFIU9DtcOl9Tid4fkp1ltndx9FJ85jscByaRYX0mw260DV2xvudZxwsre4lrbOXnInjmXGpBSzw4l4kqCi0JQJyWSkJXChvYfyqvD7FDtYQVkD/Q4neVPHkZIUZ3Y44hI2LTUS1M7Cavr6nSZHMzrb9p4C4LpVU2VroyCQBBWFLBbLQCm2p/ItnO0/5p7emyfTe+Fg+qQUpmYl09rROzA1Gw5qGtopKj9PbIyNDRHesiZUSIKKUp4EdSDME5TD6RpIsislQYUFi8XCtaumAvDSW6dMjWU0tu09DcBViyaRlBBjcjTRQRJUlMqfmU6s3Ur52Qs0t3WbHc5l06ebaO3oZeL4MUzOlOaE4WLT0hxi7VYKyhrC4pqovn4Hrx0wqveuWz3V5GiihySoKBUfax9oR3GoJHwv2vWe3pM1gfCRnBjLWnefqJf3nTY5mkt7q7iW1g6jOGK29H0KGklQUWxgHSqMy833DSQoqd4LN9etzgXg1f1n6HeEdrHEtreMJHrd6lz5IBREQU1QSqlrlFKHlFKtSqlypdS97ttjlVK/U0o1KaUalFJfCWZc0cqToAp0fci/QQylpqGdqvp2xiTEMHea7B4RbubkjiNnQjLNbT0DI+FQVFXfRnHFeeJibWyUtu5BFbQEpZTKAZ4AvgOkAh8EHlBKvQv4JqCAGcBy4G6l1F3Bii1aZbnXbTq7+yk5GX67m+89arypLc3LxG6TyYBwY7FYuM5TLLE3dKf5nt99EjCKI8ZIcURQBfO3Ohf4m9b6Ka21U2t9ANgOrAXuBr6rtW7WWp8CfgTcG8TYolY4l5vvcXcz9bQUF+Fn07IcYuxWCsrqQ7JYor2zl1f3G8URN62fbnI00SdoCUprvUtr/UnP90qpccB6oACYCBz3OrwUWBCs2KKZZ+eFfSE8xTKU+uZO9Jlm4mJt0j03jHmKJVwueGHPKbPDeYcX3zpFd6+DRbMzmJYtO0cEmynzIkqpFOBZYB9wyH1zp9chnUBisOOKRnOnjSc5MYbqhnbO1LWaHY7P3iquBWBZ3gTiY+0mRyOuxM3ukcm2t07R0dVnbjBe+vqdPPemMb13y4YZJkcTnYKeoJRSs4G9wDngVqDN/aMEr8MSgfYghxaV7DYrK+cZLRA8b/rhYPcRY3pvrUzvhb1ZOWksmJFOV09/SK1F7Sqspqm1mylZyQMbLIvgCnYV31UYo6angVu11t1a62agDqNIwiOPi6f8RACtzjcS1J6i8EhQjS1dlJxqItZuZak0J4wI791ojFCe3VUREvvzuVwuntlRAcAtV82Q0nKTBG1uRCk1A3gO+KrW+heDfvww8HWlVBGQBNwP/DxYsUW7RbMySIizU1nTQl1jB1njx5gd0og8I70leZkkxktVVSRYmjeBnAnJnD3Xxq7CKjYvm2JqPEXl56msaSE1KY4NUlpummCOoO4DkjFKy9u9vr4PfA04ChwDDmCUoz8YxNiiWmyMjeXuar5wGEXtLpLpvUhjtVrY6h5FPbW9wvRmhk++UQ7ADeumERsjDTDNErQRlNb6i8AXRzjkPveXMMGa/Gx2FlbzVnENWzfNNDucYTW3dXOsshG7zcpyaU4YUTYsmczDL5ZwqraVAt3AEpOqM49VNnJY15MQZ+N6924XwhxydaMAjOmyWLuV0tPNNLZ0mR3OsPYW1+JywWKVIRdNRpgYu40b1xkVff94rcyUUZTL5eLhF0sAeM9VM6W/mMkkQQkAEuLsA59Y94ZwNd+b7uq9NQtkei8SXb9mGkkJMRyrbORQafA3MS7QDRyrbCQ5MUZKy0OAJCgxwLMjw54QTVD1zZ0UV5wnxm5l9YKJZocjAiApIYbbtswG4M/PHcPhDN4oyhg9GcXD79s0S0boIUASlBiwfG4WdpuFoxXnudDWY3Y477DjcBUul9GYUN48ItcNa6eRkZbA6bo23jh4NmjP+1ZxLeVVLaQlx3HDumlBe14xPElQYkBSQgyLVSZOF+wsqDI7nIu4XC5ed79ZbV4m7bYjWWyMjQ9fNweAR7aV0NPnCPhzOhxO/rqtFIDbt8yW3UlChCQocRHPm//rh4L3ydUXJ85eoKq+ndSkOBbLVf0Rb+OSyUzLHsv5lm6e21UZ8Od7dlclZ8+1kTkukWtX5Qb8+YRvJEGJi6yYa0yfVVS1cLo2dPbme8OdMK9aPElaa0QBq9XCPTfMA+Dx18pobu0O2HPVnu8YGD19ams+MXZ5fYUK+Z8QF4mNsbF+0STg7aRgtr5+JzsLqgGjPYOIDotVBkvyMuno7ucXjxcGpOzc5XLxq38W0tvnYOOSyQPtZ0RokAQl3mHzUiMJvHGoKqhVVMM5XHqO1o5eciYkM2OStDyIFhaLhc+8fxFj4u0cOH6OV9x9mfzptQNnOXLiPMmJsXz8PfP9/vjiykiCEu+Ql5vGxPFjaGrt5siJBrPD4Y1DRsHG5mU5smlnlElPTeCTW/MB+MMzxX5tatjc2s0fnz0KwCdumS8X5YYgSVDiHSwWy8BUWjDLfIfS0t7DvmN1WCzGwrmIPhuWTGbtwmy6ehz87O8FfhnV9/Y5eOChA7R39bEkL1M2hA1RkqDEkDYtNX5h9xTX0tltXhO5l/edpt/hZInKJD014dJ3EBHHYrHwqa35pCXHcayykT8/d+yK1qOcThc/+3sBJaeaSE+J57O3LZKReYiSBCWGlDV+DPOmj6e3z8GuwhpTYnA4nANtwD17tInolJIUxxc+uASb1cLTOyr420v6sh/rkZdK2VVYTUKcna99fBXjU+SDT6iSBCWGde1KoyfPc29WmrJx575jdZy/0EV2+hjpaCpYrDL50oeXYbXA31/RPP5a2agf4/ndJ/nHq2VYrRa+fNdypmVL0U0okwQlhrV+0SRSk+M4VdtKccX5oD//c2+eBIytb6xWmYIRsHZhNl/44BIsFvjLCyU8sq3Upw68XT39/PTRwzz4ZBEAn9yab1o7D+E7SVBiWDF2G+9eY+xJ9syOwF/N782TFONjbVy93NzuqiK0bFyaw6ffvwgwRlKf/fEbHCkbvtq0ouoCX/jpDl4/eJbYGBufuW2R9HkKE7LhlBjR9atzefy1Mg6U1FHT0E52RlJQnvf53cboafOyHNkYVrzDtSunkpGawINPFlFV385//3YPS/IymZ2TxuTMJFKSYimuaOTA8TpO1hg7okzNSuY/7lzGlKyxJkcvfCUJSowoNTmOjUsm88r+M/xrVyX3uq9JCaT2zt6BXSykOEIMZ7HK5Jdf2sRT2yt47NUyDpfWc3iIHlLxsTa2LJ/CPTfNI07at4cVSVDikm6+agav7D/DqwfOcMf1c0gK8IjmhT2n6Ol1sGhWBjkTkgP6XCK8xdht3LZlNlcvz6GwrIHqhnaq6ttpaulm1pRUls/JYv6M8cRKYgpLkqDEJeVOHMvCWekcOXGel/eeZuummQF7rvbOXp7cXg7ArZtnBex5RGQZn5Iga5URSIokhE9uvspof/3MznK6evoD9jxPbi+no6uPhbPSWTg7I2DPI4QIfZKghE+W5U1gVk4qTa09PO0e4fhbc2s3z7p7/9x5/ZyAPIcQInxIghI+sVotfOxmY7fnJ7aX09jS5ffn+MdrZfT0Olg5Lws1dZzfH18IEV4kQQmfzZs+ntULJtLT6+ARd4M3f6lv6mTbW6ewWODDMnoSQiAJSozSPTfMxWa18OqBM5ysafHb4z78Ygn9DhcbFk8md6JcpyKEkAQlRik7I4kb1k7D5YL/+9eV7Srtsbuohu2Hq4i1W/nQu/L8EKUQIhJIghKjdvs1ijEJMRSWNfDS3tNX9FiNLV386vFCAD5y0zwmpo/xR4hCiAggCUqM2tgxsdz73gUA/PapYsrONF/W43j68rR1Gk3jblg7zZ9hCiHCnCQocVk2Lc3h3Wty6Xc4eeChA7S094z6MZ7bXUlhWQPJibF87vbF0jROCHERSVDisn38PQtQU9M4f6GLH/310Khace8uquFP/zoOwGduW8i4sfGBClMIEaYkQYnLFmO38uW7lpOSFEvhiQa+/5cDtHdduj3887tP8v2/HKDf4eSWDTNYvSA7CNEKIcKNJChxRdJTE/jyXctJjLfzVnEtn/vJ9mHXpJxOF399sYQHnyzC5TJ2i/joTfOCHLEQIlzIZrHiis2fkc7PvrCRHzx8gPKqFv7zl7u4ZuVU1JQ0Zk5OxeF0seNwFTsKqmhs6cZqgfvev4hrV041O3QhRAiTBCX8YmL6GH7wmfX86bnj/GtXJS/uOcWLe06947jMcYnc+94FrJibFfwghRBhRRKU8JsYu41P3LKAqxZPorj8POVVF6ioaqG3z8GqBRPZuGQyc3LHSbWeEMInkqCE3+VNHUeebPYqhLhCUiQhhBAiJEmCEkIIEZIkQQkhhAhJkqCEEEKEJElQQgghQpIkKCGEECFJEpQQQoiQFAnXQdkA6urqzI5DCCHEKHm9d9sG/ywSEtREgDvuuMPsOIQQQly+iUCF9w2RkKAOAOuBWsBhcixCCCFGx4aRnA4M/oHF5fK9yZwQQggRLFIkIYQQIiRJghJCCBGSJEEJIYQISZKghBBChCRJUEIIIUKSJCghhBAhSRKUEEKIkCQJSgghREiKhJ0kfKKUWgg8COQDlcBHtdbvuHJZKTUF+COwCqgHPqO1fiGIcV4DfA+Y5X7+H2qtfzvEcZuBV4Aur5u/r7X+dpDi/CjwW6DH6+b7tNYPDTrOtPOplLrDHaO3BOA1rfW1g4415XwqpVYAz2mtM93fxwK/BG7F2BnlJ1rrB4a5rwX4NvAJIBb4E/AlrXV/EOLMBH4OXA1YgBeBz2mtm4e5/xlgPODZGaBaa62CEGcc0Ab0eh22Z/D/v/tYM89n+6BD7EAcMElrXTPE/QN2Pod7DzLjtRkVCcp9Yp8BfgZcBbwPeFkpNVVr3Tro8L8DbwE3AOuAp5VSi7TWlUGIMwd4ArjbHe9S4CWl1Cmt9UuDDl8CPK61/kCg4xrGEuDHWusvX+I4086n1voR4BHP90qpxcDLwJeGODyo59P9C/wx4EeDfvRNQAEzgBRgm1KqWmv9lyEe5hPAVozYe4CngP8CvhWEOP8AtADTgBjgYeBXwIeGeIx0YBIwVmvd4a/YfIxzAdCktc7y4WFMO59a6ySvY+zAG8D2YZJTwM7nSO9BwEaC/NqMlim+jUCM1vpnWus+rfXfgWPA7d4HKaVmA8uAr2mte7XWrwPPYryggiEX+JvW+imttdM9wtsOrB3i2KVAYZDiGsolnz8Ezqd3LDEYyeobWusjQxwS7PP5TeBTwHcG3X438F2tdbPW+hTGG9m9wzzG3cDPtNZVWusG4BsjHOu3OJVSVsAJfFNr3aG1vgD8HuMDyFCWAicClZyGi9PruX39fzXlfA7hPzGS/teH+Xkgz2cuw78HBf21GRUjKGAuUDLotlKMT1eDjzsz6D++FFgRwNgGaK13Abs83yulxmFshPvwEIcvATKUUp/CmGJ5DPhvrXXPEMf6lVLKhjFVeqdS6idAJ8Yn6u9rrb03dzT1fA5yH8b03a+H+Xmwz+eDWuuvKaU2em5QSqVibJp53Ou4oV6nHnOHODZbKTVOa90UqDi11k7glkHH3QIUDPMYSwCrUmo/xojrMPB5rfXg30m/xun13JlKqSJgArDT/dzVQzyGKefTm1IqG2OksdZ9nocSsPN5ifegoL82o2UElYTxJuqtE0i8zOMCTimVgjHa2Icx1Pb+mR2owhg2zwE2A1sw5nyDIQM4CDyE8QtyK8anwk8NOi4kzqd7ivdLGKOnd+yObMb5HGrqBuN8wcXnbKTzNfj8ev7ut/M7TJwXUUrdj5Gg/nOYQxzAfowpn6kYiexFpVQw4uwAdmOslSmMDylPDXNsKJzPLwDbtNYjjfoCfj7hHe9Bh9w3B/W1GS0jqA6MxXFvicDghUlfjwso99TYMxifQO4Y/EnKvdB4tddN5Uqp7wLfB/4j0PFpreuADV43FSqlfoGxtuc9QgmJ8wlchzEl9fxQPzT7fHrxjDS9z9lI52vw+fX88gfl/LqnTX8B3ARs1lqXDnWc1voHg+73FeD/YUxV7RrqPv6itf7ioOf+ItCglMrRWp8ddLjZ59OGMTU2YnO7YJzPwe9BvH1egvrajJYR1HGMT0/e8rh4COo5bopSKuESxwWMUuoqjE8sTwO3aq27hzhmklLqR+6RgUcs8I5jAxTjPKXUNwfdPNTzm34+3d4D/GO4KROzz6eHuwKujotfqyOdr8Gv6zyg1r0mFFBKqWSMqsflwIqRPvErpT6vlPJen7JhfDgO+PlVSn1LKTXH6ybP//FQz23a+XRb4/7ztZEOCvT5HOo9yKzXZrSMoN4ALEqpL2CUSb4PYw3loqG+1lorpY4A33V/KlmD8ea2OhhBKqVmAM8BX9Va/2KEQxsxPtV0KqW+hTHN9t/A/wU+SgAuAP+ulKrCKCFfDHwW+LT3QWafTy+rgP8Z4edmn09vDwNfd6+ZJAH3Y5RzD3fs/Uqp1zA+sX6DodcrA+HvGB9w12utB0/jDpaLsV55I8Zr5/vACYy1k0DLB5YppTzVhT8Hnncv3A9m5vkE43W6d4S1J49cAnQ+L/EeFPTXZlSMoLTWvcD1GImpCfgqcIvWukEpdcegaxDeh7EOUY+x8P8xrfXRIIV6H5AMPKCUavf6+r53nO5R1fUYJfONGAu/jwM/CUaQ7gXmmzGqcloxylK/rbX+Z4idT49c4KK5/1A6n4N8DTiKUWV6AOPcPuj5ofv14JkCehAjzj0Yb1DH3fcPKKVUPvBujGKXeq/XadUwcX4Z2IuxVlIPTAdu0loHowP2x4BmoBw4hXE91J3DxGnK+fSSy6DXqUcQz+ew70GY8NqUjrpCCCFCUlSMoIQQQoQfSVBCCCFCkiQoIYQQIUkSlBBCiJAkCUoIIURIkgQlhBAiJEXLhbpC+J1S6s8YW9MM55sYO0G/ASRrrYO5Zc5u4C6tddkIx1kxrqe5U2utgxGbEKMhIyghLt/nMHZ4nojR0gWMi1c9t/0I40LFiby9z14wfBY4MlJygoFdyb+F18WWQoQSuVBXCD9QSs0HioFp7l45ZsURD5zB2LzVpx07lFIVGDt8bA9kbEKMlkzxCRFA7r4/A1N8SikX8EHgKxibaR4EPozRDuROjK2jvqK1fth9/2TgxxgtTVzA6xit1Ydr2/AB4IJ3clJK/Q9Gh9MMjL5o/6W1ftHrPk9hjAa3++GfLITfyBSfEMH3PeDzGJuDTsHY5LMVY2fwJ4HfKqU8vaF+h5HI3oXR4sSF0YJ7uA+XNwDbPN8opd7rfq4PY+wo/TzwuFJqrNd9tgFbRnhMIUwhCUqI4PuV1voNd4uK5zB65PyXu1DhJxh9dKYppaZjjIg+pLU+4B4V3Ymxqeh1wzz2MozNPD1ygR7gtHvq8VsYje76vI45jrE7dZ5f/nVC+Il8YhIi+Mq9/t4JnPLq9Ovp6ROH0S0VQCt1UTuzRIxR1XNDPPYE4LzXhCQfmAAAAXVJREFU93/FqDSsVEodwuiQ+ietdZfXMY3uPzNH+e8QIqBkBCVE8PUN+n64/j9297GLgUVeX7OBPw1zHydg8Xzj7nu0FGPEtQe4ByhyF3V4eN4HgtH+QgifSYISInSVADHAGK11uda6HKgFfoiRpIZSh1EMAYBSaitwr9b6Za315zBGXm0Y/Zw8MrzuK0TIkCk+IUKUuyPxs8BflFL3AQ3AdzGKK0qHudshYKHX9zbgh0qpcxgVg6uALPffPRbydlM/IUKGjKCECG13YySTpzG6mKYA12itLwxz/PMY1X4AaK0fB76OMeoqA74DfFpr/brXfa4CtgWpw60QPpMLdYWIIEqpRIzW5tdprQ/7cLwVOI1RKbgrwOEJMSoyghIigmitOzFGS/f5eJf3AJWSnEQokgQlROT5KZCvBtWmD+YePX0V+GRQohJilGSKTwghREiSEZQQQoiQJAlKCCFESJIEJYQQIiRJghJCCBGSJEEJIYQISf8/vx6/Sfj8upwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_position(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After reaching the lowest point, the jumper springs back almost to almost 70 m, and oscillates several times.  That looks like more osciallation that we expect from an actual jump, which suggests that there some dissipation of energy in the real world that is not captured in our model.  To improve the model, that might be a good thing to investigate.\n",
    "\n",
    "But since we are primarily interested in the initial descent, the model might be good enough for now.\n",
    "\n",
    "We can use `min` to find the lowest point:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "5.413731910623594 meter"
      ],
      "text/latex": [
       "$5.413731910623594\\ \\mathrm{meter}$"
      ],
      "text/plain": [
       "5.413731910623594 <Unit('meter')>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "min(results.y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the lowest point, the jumper is still too high, so we'll need to increase `L` or decrease `k`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's velocity as a function of time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_velocity(results):\n",
    "    plot(results.v, color='C1', label='v')\n",
    "        \n",
    "    decorate(xlabel='Time (s)',\n",
    "             ylabel='Velocity (m/s)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdeZxbZ3no8d/ROtJoFm/jfV9eb+M9JnESxwkJSSghBFIgDSkU2kspt6Ut0ALlsqcUmtJ0ZSlctlsaCGFLQvbYZHUSO47H6xvv9ngZ22PPImmk0XLuH0eS5cnMWB5LOkfS8/18/IlHc3TOkxNFz3m35zVM00QIIYRwGpfdAQghhBCDkQQlhBDCkSRBCSGEcCRJUEIIIRzJY3cApaaU8gOXAceBlM3hCCGEOJ8bmAi8orWO5/+i6hMUVnJ61u4ghBBCDOtq4Ln8F2ohQR0H+O///m8mTJhgdyxCCCHynDhxgjvvvBMy39X5aiFBpQAmTJjAlClT7I5FCCHE4N4wBCOTJIQQQjiSJCghhBCOJAlKCCGEI0mCEkII4Ui1MElCCCEqRjqdpr29nUgkYncoReH1emlpaaGxsfGi3ysJSgghHOT06dMYhoFSCpersju5TNOkr6+Po0ePAlx0kqrsf3shRshMJYgd2U2iqwPZckY4SVdXF+PHj6/45ARgGAbBYJDJkydz8uTJi36/tKBETUnHIvRseYLulx8mFT4DgOEL4GuZTsOSa2lY9mYMw7A5SlHLUqkUXq/X7jCKKhAIkEgkLvp9kqBEzeh59XE6n/oRZn8fAJ7m8ZiJOKlIF/H23cTbdxNr14y9+U9weXw2RytqWbU9JI3030cSlKgJvds2cPqRbwNQN30xzZffSmD2cgzDIBnuIrpnE52Pf49w29MkTh9h/Ls+iadxjM1RC1HbKr+TU4gLiO59lVMP/ScAo6//AJPe90WCc1bknuo8oWYal1/PpPf/PZ6mccSP7eHo9z9FMtxlZ9hC1DzHJCil1A1Kqc1KqR6l1F6l1Iczr/uUUt9RSp1RSp1SSn3a7lhF5YgdfZ2OX9wD6RRNV7yD5jfdMuSx/gkzmfzBr+OfrEiFz3DqwX/HNNNljFYIkc8RCUopNRV4APgK0AzcAXxVKXUj8EVAAbOxts54v1LqD+2KVVSOVF+Yjvv/ATMRJ7RkHaOvfd8F3+MONjL+nR/HFQjRt38LPa/8tgyRClF5Pv7xj3P33Xfnfk6lUqxZs4aXXnqpaNdwyhjUDOAnWutfZn5+RSm1AbgSeD/wAa31WeCsUuoe4MPAj+wIVFSOs8/+jFSkm7qpCxj31o8UPFDraRzDuN/7Mzp+/nU6n/4xddMX4x8/o7TBCjGE4/fdTd++V8t2vcDsFUx8799d8LjbbruNT33qU3zqU5/C7Xbz/PPPU1dXx+rVq4sWiyNaUFrrZ7XWf5r9WSk1Gmvzqi1YOy3uzDt8N9Ba3ghFpek/3U7P5kfBcDHmxj/GcF/cs1i9ehMNy98CqSQnf/XPpBPxC79JiBqyZs0aXC5XrsX04IMPcssttxR1BqJTWlA5Sqkm4DfAS8DmzMvRvEOiQLDccYnK0vnEDyCdomH5DSNu/Yy54QPEDu8gcbqd7pceZNRVtxc1RiEKUUhrxg4ul4u3v/3tPPjggyxfvpwnn3ySX/ziF8W9RlHPdomUUvOAjUAHcDvQm/lVIO+wIBAuc2iigkT3bqZv/xZc/iCjr7ljxOdxef2MvelPAOje+GtSffKxEyLfbbfdxpNPPsnjjz/OvHnzmDlzZlHP75gEpZRai9Vq+hVwu9Y6lhl3OoE1SSJrPud3+QmRY6aSdD75AwCar/593PVNl3S+wIxW6qYvJh2P0v3Sg0WIUIjqMXv2bKZPn869997LrbfeWvTzOyJBKaVmAw8Bn9Naf1prnV8c7cfA55VSY5VSM4BPZF4T4g3C258h0XkM7+hJNK26uSjnzLbCul95iFS0pyjnFKJa3HbbbZw6dYq3vvWtRT+3IxIU8FGgAWtqeTjvz9eAzwHbgR3AK1jT0b9lX6jCqUzTpDszLbx5zW0Y7uLUM6ubOp/A7OWY/TG6XvxVUc4pRLW488472b59O83NzUU/tyMmSWit/xr462EO+WjmjxBDirdr+jsO4Ao2Ur/oqqKee/Q1d3B03xZ6Nj1C0+pb8DSMKur5hRBv5JQWlBCXrHuT1XpqXHZ90Yu9+ifOJjhvNWayn+6N0ooSohwkQYmqkOw9Q2T3RjBcNK68sSTXGHXV7wPQu/Vp0v2xklxDCHGOJChRFXpefQzSKerVajyNY0tyDf/EWfgnzSUdjxLe+XxJriGEOMcRY1BCXAozmaB3yxMANK4q/kyifI0rb+TUsT30vvoYjcveXNJricqWDJ+ld+vTRPe+iruuHk9zC56mcQTnrsI3ZvKw7zVNs6r2hBrprtWSoETFC+9+kVSkG1/LNOqmLSzpteoXrKHziR8QP76P+LG9+CfNKen1ROWJdxyk67n7ibz+CqRTb/j9mfU/oXnNO2i+8l2DjpW63W4SiQQ+X/VsmtnX1zeiXYIlQYmKF25bD0DjyptL/tTp8vppWLKO7pcfoufVxxgnCUrkiR7YSsf9X8dMxMBwEZy3moYl68A0SXSfov/4PsI7nqXruZ8T2fUCY9/6EQIDHqqam5vp6Ohg8uTJuFyVPQpjmiZ9fX0cPXqU8ePHX/T7JUGJipYMd9F3cDu4PNQvWFOWazasuJHulx8ivOM5Rl//Adx19WW5rnC2yO6NdPzqnyGVpH7RVYx58/vxNIx+w3GNK2/k1MPfJNF5lOP//QUmvPvTBGcvz/1+7NixtLe3o7UuZ/gl4/V6GT9+PI2NjRf9XklQoqJFdr8IZprg7OW4A6GyXNM3ZhKBGa30HdxGeNsGmi77vbJcVzhX79anOfXwN8FM03jZWxlzwx9hGIO3fuqmLmDKH/8TnU/+gJ7Nj9LxwD1MuuvL+CfOAqwirNOmTStn+I5V2e1HUfMimdl0xV6YeyENK6yp7D2vPj7iAWBRHWLtu3PJadTa9zDmhg8OmZyyDI+XMTd+iNDitZiJGCd+ejeJro4yRVw5JEGJipXs6SR2ZBeGx0f93MvKeu36eZfhCjaSON1O/4kDZb22cI50fx8nf/2vYKZpuvxWRl397oLHQQ3Dxbi3/RmBGa2kIl2c+J+vkOrrvfAba4gkKFGxwrus1lNwzkpc/sAFji4uw+0hlBnzCu94pqzXFs7R+eQPSXZ14GuZweh1F7+1i+H2Mv72v8HXMoPEmWN0PvnDEkRZuSRBiYoV2ZHt3rvSluuHFq8FILzjOcxBphOL6hbZs8laf+f20HLrx0ZcnNjlD9Lyzo9juL2E29bTd3BbkSOtXJKgREVKnD1B/PheDF8dwdkrbInBP3kenubxpMJniR3aYUsMwh6paC+nH/4mAKPX/QG+lkub1OAbM4nmzI7Np377LdKJ+CXHWA0kQYmKlC01VD9vNS6v35YYDMMgtOhqK54dz9oSg7BH1wu/IBXpom7qAppWv60o52y+4la8Y6eQPHuCrucfKMo5K50kKFGRIjtfACC0sLyz9wYKLc4kqN0bSSf7bY1FlEey9ww9mx8FsGbsudxFOa/h9jLu9z4CQNeLv6L/5OGinLeSSYISFSfZfYr+kwcxvHUEZi6xNRbf2Cn4xs/EjEeJ7t1sayyiPLqefwAz2U/9/Mtza5eKpW7KfBpX3AjpFJ1P/aCo565EkqBExYnufRWAwMwlGJ7i7Jp7KXKTJbZLN1+1S3SfpGfLk4DBqLXvKck1Rq27A8MXoG//VmLt1VFNYqQkQYmKE91nJajgnJU2R2IJLbwSMIju3UyqL2x3OKKEup69H9JJQouvxjeuNNUe3IEGmi6zqvKfffZnJblGpZAEJSpKOtmfm4YbnGPP7L2BPI1jqJuxGFJJInqj3eGIEunvPEZv2wYwXIy6+t0lvVbT6lsyrajXaroVJQlKVJTYoR2YiTi+8TMHLcRpl3Oz+Z6zORJRKt0v/hLMNA1Lr8M7emJJr+UOSisKJEGJCpOdiOCU7r2s+vmXg9tD7OB2kj2ddocjiiwV7c09fDRf8Y6yXFNaUZKgRAUxTfNcgprrrATlrqvPJE0zV4JJVI/etqcxk/0EZi0veespS1pRkqBEBUl0HiXZdRJXsBH/xNl2h/MGuTVR26Wbr5qYZpqezY8B0LjqprJeu2n1LRheP337X6P/dHtZr+0EkqBExchOLw/OWla0xZHFFJyzEsMfpP/EPvo7j9odjiiSvn1bSHZ14GlqOW9jwXJwBxsILb4GgJ5Nj5T12k4gCUpUDKeOP2W5PD7q1eWArImqJt2brKoRjStvtOXBqCnTauvdtoF0PFr269tJEpSoCOl4lNiRXWC4CMxaZnc4QwottkovhXc8KxsZVoHE2RP07duC4fbSsPQ6W2LwtUynbtoizP4Yvdt+Z0sMdpEEJSpC36EdkE7hnzy3bFu7j0Rg+mLc9c0kz54gfmyv3eGIS9Tz6mOASf2iK3EHG22Lo3HVzVY8mx6pqQcfSVCiImQX5wZmtNocyfAMlzu3/Xx4u2xkWMnMVNJamAs0rijv5IiB6uddhjs0mkTnUWI1tF+UJChREWKHsgnK3uKwhWjI1OaL7HpeNjKsYH37t5KO9uAdNxX/pDm2xmK4PTSueAsA3TU0WUISlHC8VKSb/pOHMTw+6ibPszucC/JNmIV3zCRSkW76DrTZHY4Yod7t1nhPaNFaDMOwORpoWH49uDxE92wi2XPa7nDKQhKUcLy+Q9sBqJu6wBHVyy/EMIzc1GDp5qtM6Xgf0ddfAc5NfLGbJzSKerUazHTNTJaQBCUcL9sKCcxYbHMkhQtlxqEi+iXS/X02RyMuVkS/hJnsp27qArxNLXaHk9OwZB0A4bYNNTFZQhKUcLxcC6oCxp+yvKMm4J+iMBNxIpkncVE5si3f7F5fThGYtQx3fTOJM8eIH9tjdzglJwlKOFqi+yTJsydw+YP4J8y0O5yLkp0sEd4m3XyVJNl71po16vZQv2CN3eGcx3C5c0mzt229zdGUniQo4Wixg5nW0/RFjixvNJz6BVeCy03fga0kw112hyMKFN75HJhpgrNXOHLNXbabL7LzedLJfnuDKTHHJSil1Gql1Mm8n31Kqe8opc4opU4ppT5tZ3yivCpl/dNg3MEGq3abmSayUwrIVopc916rs7r3snwt0/GNn0k6FslN5KhWjklQSilDKfXHwOOAL+9XXwQUMBu4DHi/UuoPbQhRlJlpmhWdoODcGIbM5qsM/Z3H6D+xH8MfdGzNRzjXisouJK5WjklQWInoI8BXBrz+fuBurfVZrfVB4B7gw2WOTdgg0XmUVPgs7vpmvGOn2h3OiATnrsLwB4kf31eT2yVUmqjeCFiVG1we3wWOtk9o0dVW9/H+10iGz9odTsk4KUF9S2u9EtiUfUEp1QxMBHbmHbcbqMzHaXFR+vLHnxywUHIkXF4/ofnZCufSinK6iH4ZgHr1JpsjGZ67vong7BVgpgnvqN7K+Y5JUFrrY4O8nB2hzK8xHwWCpY9I2C12xHouCUxbZHMklya/m8800zZHI4aS7OkkfmwPhsfn6Ir5WaElmcXgVTxL1DEJagiRzD8Dea8FgbANsYgyMk2T2OFdANRNW2BzNJembvoi3I1jSXafInZkt93hiCFE9EsABGYvx+X12xzNhdXPWYWrrp7+jgP0nzxsdzgl4egEpbU+C5zAmiSRNZ/zu/xEFUp2dZAKn8EVCOEdO8XucC6JYbhoyG4HXyMlaipRNkE5vXsvy/B4c+u0snUDq42jE1TGj4HPK6XGKqVmAJ/IvCaqWOxIpvU0dQGGUQkf0+Fla/NFdr1Q9WtXKlEq2kPs8E5wuR09e2+ghtZzNR+rsXJ+Jfyf/zlgO7ADeAV4APiWrRGJkst1701daHMkxeEbN9VauxKP5rauF84R3bMJzDSBGYsduTh3KP4p8/E0t5DqPUPs0A67wyk6j90BDKS13gA05/0cAz6a+SNqRH4LqlqEWq/hTMcBwtt+R2j+FXaHI/LkuvfmVUb3XpZVOX8tXc/9nN7tvyMws3LqVRaiElpQosYkw10kzhzD8Porrv7ecEKLrgLDRXTvFlLRXrvDERnpeB99+7cCBsF5q+0O56Jlu/kiuzeSTsRtjqa4JEEJx4m1Z1pPk+dhuB3XyB8xT2gUgZlLIZ0kvPN5u8MRGdH9WzBTCfxT5uFpGGV3OBfNO3oS/klzMftjRDPruKqFJCjhONU2/pSvIbd2ZYO9gYic6B5rTLC+AltPWaFMK6raNjIs6PFUKdUK3AysAlqAFNb071eAh7TWe0sWoag5ufGnCl//NJjgvNUYvgDxY3vo7zyKb8xku0OqaaaZJrrvVYCKmr03UGjhVXQ+8YNM5fyzeEKV1xIczLAtKKXUWqXUeqzyQ7cAZ4EXMj/HgPcBO5VSTyilnFn6V1SUdDxKf8dBcLnxT55ndzhF5/L6qc9MkKjmCgCVIn5sL+loD56mlopeb+cONhCcU32lj4ZsQSml/i+wCPh34Dat9aAb2iilGoE7gHuVUm1a6w+UIlBRG2LtGsw0/olzK2I1/0g0tK4l3PY04e2/Y9Q176mKdV6VKjvlPzhnRcXWe8xqaF1H9PWXCbf9juY3vd3ucIpiuC6+h7XWH7zQCbTWPcC3gW8rpW4vWmSiJlVz917WeaWPDu8iML2yaw1Wsujeyu/eywrOWYGrLkT/yYPEOw7iHz/D7pAu2ZCPblrrBy72ZFrrn19aOKLWxdqtWnXVOEEiyyp9lN0OvroGtStJsvestfeTx0ddFTwkGB4voYVXAhCuktJHBfUtKKWCSqkvKqXmZn7+tlIqrJR6Wik1sbQhilphppLEj1nzbeqmqAscXdmys67Cu1+surUrlSK6z+reC8xorZru5FBmI8PwtuoofVRo5/e/YE2I8CmlbsXaRPCvgD7gX0sUm6gx/ScPYSbieEdPxB1stDuckvKNnYJ/4hzMeLTqt+12qmrq3svyT5qLd/REUpEu+g602R3OJSs0Qd0K3KG13gH8PvCE1vq/gE8CbylVcKK2xNo1AP4qbz1lnVu7ssHeQGqQmUzQd2ArYI3dVAur9NG5ArKVrtAEFQA6lFIu4Ebg0czrJtaaKCEuWW78aXKNJKhFV2W27d5a1dt2O1HfkZ2Y/TF8LdPwNI2zO5yiCrVa45uR3RtJx/tsjubSFJqgXgH+Fvg8MAr4pVJqEvBlYGOJYhM1Jn70dQDqpsy3OZLycAcbq3LtSiWoxu69LG/zeOqmLcRM9hPZ/aLd4VySQhPU/wbWAB8D/iyzPfunsTYS/IsSxSZqSLKnk2T3KQx/EO+4yl0webEaWq8FINxWHbOuKkVfpnpEYHb1dO/lq5bu4+EW6q4BNmqt01rrncCyAYd8RmstJZlFUcSyrafJc2tq4Wpwzgpcgepau+J0ie6TJDqPYfiD1FVhtRKA0II1dD72PWKHdpDoOom3ucXukEZkuG+Ce4BTSqlfKaX+TCk1J/+XkpxEMcVz40+10b2XZa1duQqQArLl0rfvNcCaXl5N1fLzufxBgsoqflvJkyWGW6i7BpgJ/BBYDDyqlDqglPqOUup2pVR1VCMUjlBrM/jyhVrXARDe/mxVrF1xuuh+K0EFZw3sFKouDdnP1bYNmKZpbzAjNOzjQ6aM0S8zf1BKzcSaVv5e4FtKqX1YU84/W+pARfVKJ/uJnzgAGNRNnmt3OGXnnzQH75hJJDqP0bd/a1VNe3YaM50idnAbAIEqT1CBmUtwh0aROHOc+NHXK3Lx+0V19mutD2itv621vh1r242/wKpqLsSI9R/fD+kkvpapuPxBu8MpO8Mwcq2oSh/Udrr40T2k41G8oydV7LhMoQyXm9Diq4HK/VwV3AGrlFqHVd18YE2QaDEDErUnu/7JX2PjT/kaFq/l7Ib/IapfJhWL4K6rtzukqhTdvwWo/tZTVkPrOro3/obIzucZc8Mf4fL47A7pohS6YeG9wJ8Dh3lji8kEvlHkuEQNyY4/VWIXRLF4msZRN2MxsYPbiOx6gcblN9gdUlXq25+pHlEjCcrXMh3f+Jn0dxwgumcToQVr7A7pohTagvpD4INa6x+WMhhRe0zTJH40m6Cqc8pvoRpa1xE7uI3etg2SoEogFe21ihG7PVVRvbxQDUvW0fnEAcJtGyouQRU6BhUFXi5lIKI2JbtPkop04wo24hlV24Xx6+e/CcNbR7x9N4kzx+0Op+pYtfdM6qYuwOWrszucsgktuhoMF9F9W0iGB9131rEKTVBfAe7JzOITomhy5Y0mza34HU0vlcsXoH7+5UDlDmo7Wa1MLx/IXd9EcPbyiiypVWiC2gVcAexVSqUG/ilhfKLKZStI+Kt0Rf/Fasjfz8dM2xtMFTFNMzf+VCsTJPKFlmRKalXYBpmFjkF9B6so7PeRWXuiiOJH9wBUbcmZi1U3fRGexrEku0/KdvBFlDh9hFT4DO76Znwt0+0Op+zq566ytoPvOFBRJbUKTVBTgZu11vtLGYyoLfkLdP2T5lzw+FpgGC5CrdfQ9fwD9LZtkARVJNnN+wIzl9RkV3J2O/ieVx8jvG0D/vEfsDukghTaxfcEsLaUgYja03/CWqDrHVebC3SHkl20G9n9Aul+WQdfDNFs997MpTZHYp/cdvAVVFKr0BbUS8C/K6XeBewFEvm/1Fr/TbEDE9XvXAVz6d7L5xszCf/kecSPvk5Ev0RDZusEMTJmKkHs8E7AakHVKv+kuXkltV6riL2wCm1B3YC1aWEIa9uNy/L+rCpNaKLaxXMTJGqv/t6F5Bf6FJcm1v46ZiKGd9xUPA2j7Q7HNueV1GrbYGsshSqoBaW1vrbUgYjaE5MJEkOqX3glnU98n74D20j2nMbTONbukCqWtf6ptrv3shpar7FKar3+Cqm+MO5AyO6QhjVkC0op9XmlVKDQEymlGpRSXypOWKLaJXs6SfWctnbQHVs7O+gWyh0IEZy3CjDp3Va5+/k4QXaCRLCGu/eyPI1jCcxYjJlKENn1gt3hXNBwXXzdwA6l1NeVUpcPdoBSylBKXaaU+hdgJ1BZy5SFbeLHMq2nSXNqagfdi5HbDn7b+ordz8duqb4w8eP7wOWhbprMiIRzkyUqoZtvyC4+rfW9SqmfA38DPK6USmIt2D0NGMA4rOrmBvAD4Eqt9eGSRyyqQm6B7iTp3htKYPYy3PXNJDqPET+2R7pCR6Dv0DYw0zVX3mg49epyTvv+i/hRTX/nMXxjJtkd0pCGfXTVWrdrrf8CmAjcBTwGtAOHgIeAdwNjtdZ/LslJXIy4zOC7oPz9fMIV8LTrRH37z61/EhaXr476+VcAzp+EU+gkiQjwcOaPLZRSS4FvAUuA/VjV1V+xKx4xcmYqaXW7ICWOLiTUuo7ulx4kvPN5Rt/wgYrbz8duMkFicA1L1hFuW0/vtt8x6pr3Orab3ZlRDaCU8gG/Bn4KNAN3Y3U7NtoamBiR/pOHMJP9eEZNwB1ssDscR/OPn4Fv/EzSsTDRPZvtDqeiJM6eINnVgasuhH/iLLvDcZS6aQvxNLWQ6jlN7OB2u8MZUkUkKGAd4NVa36u1Tmit7wN2AO+xNywxErnp5TW8QeHFyBWQbVtvbyAVJjt7r276IgyX2+ZonCVbUgucXTm/UhLUQqwJGvl2A62lvnD8+H4O3ftBerfLVN9iiR/LTpCQBbqFCC26GlzuitzPx07Z7r2gdO8NKvvgE9m9kXS8z95ghlBQglJKzShxHBcS4o1V1KNAyQu4paLdpCLddD3/gEz1LRKZIHFxKnk/H7uY6RR9ma6rwCxJUIPxjppA3dQFmIk4kd0v2h3OoAptQe1VSj2nlPpTpdSYkkY0uAgwcNFwEAiX+sKBGa24go0kTrfnBvbFyKWivSTOHMfw+Gpy24ORyhX6lNl8BYmfOEA6FsbT3IJ31AS7w3GsXOkjh3bzFZqgZmFNK/8IcEwp9aBS6r0XU2niEu0EBg5YzM+8XlKG20PDYquQu4wBXLpc997E2RjuQmsVi/o5mf18Th4k3nHQ7nAcT2bvFSa04AoMj4/YoR0kuk7aHc4bFJSgtNaHtdb/oLVeCqwEtgKfBjqUUj9USl1fyiCB9YChlPorpZRXKfVerOnmvyzxdYG83Sh3PIeZTFzgaDGcmBSIHRHD4yW06CpAHpQKkb//kxiaq66eoFoNOHO33ZFMkmgH9mGtRfJgVTf/kVJKK6WuKGZwWVrrfuBm4F3AGeDvgHdorU+V4noD+cfPwNcyg3QsTGTvpnJcsmplSxzJ+qeLl+2OCe94FjOVtDcYB0sn4sTadwMGgekln0dV8RryuvmcNs5eUB+LUqoeuBV4L/AWoAP4CfBZrfUOpZQL+E+sdUrTShGo1no7cFUpzl2I0JJ1nHnyB4TbNhCaX5I8XPVMM31ui3cpcXTR/JPm4B0zmUTnUaL7X6N+rux0M5jY4Z2QSuKfOFvW2RUgMHMJ7tBokmdPEG/X1E2db3dIOYW2oE4C/wGcAm7SWk/XWn9aa70DQGudBh4HeksTpv1Ci64Gw0V03xZSkW67w6lIic5jpONR3A1j8DTaMdemshmGkbcmaoOtsTjZue49GX8qhOFyE2q1xtl7HdZ9XGiC+gAwQWv9Ia31hvxfKKVaALTWv9BaV225YE+o2Zrqm07JVN8ROje9XMafRiq0+BowXET2vEKqr2qfBy/JuQkSMv5UqNwGmbteIJ2I2xtMnkIT1H3AG8oKKaWmYY1F1YTsZInerc56yqgU5yZISPfeSHkaxxCY2QqpJOEdz9sdjuMkw130nzyE4fFRN8U5XVVO5xs3Ff/EOZjxKNHXX7Y7nJwhx6CUUncAt2V+NIDvKqUGptbpWJMWakL93POn+vrHz7A7pIoiC3SLo6H1Wvr2byW8bQNNq26yOxxHiR3cBli15gyP1+ZoKktoyTrix/fS27bBGtJwgOFaUE9gLYSNZH7uy/w9+ycMvAS8o5QBOolM9R25dH8f/aeOgMuNb4IU7rwUQbUawxcgfmwP/afb7Q7HUaLSvTdioYVXgctD34E2kr3OaHcMtwiMsO0AACAASURBVGHhaeCDAEqpg8A/aq0HlhuqOaHWdfRsfpTwjmcZfd1dsti0QPHj+8BM4xs/G5fXb3c4Fc3l9RNasIberU/R27aeMdfdZXdIjmCaJn37M/X3Zi2zOZrK4w42UD9vFZHdGwlvf4bmK+xvewzXxfdW4AmtdQJ4BVin1ODVp7XWvy1NeM7jnzQH79gpJE63y1TfiyATJIqrYem19G59ivC2Zxi97g+kWjeQOH2EVPgM7vpmvONKstql6oVa1xHZvZHetvU0XX4rhmHYGs9wj/8PAROwppg/NMxxJlAz/3cYhkFD6zrOrP9/hNvWS4IqUKxdJkgUk3/KfDyjJpA8e4K+g9ukxQBEM62nwKyltn+xVqrg7OXn1R6tmzTH1niG6+JzDfZ3AaHFazmz4SdE9mwi1deLOyCLAYdjmiaxoxqQPaCKxXpQuoazz/yU3rb1kqAg170n659GLlt7tPvlhwi3rbc9QRWceJRSH1RK3Z7388+UUjXZ+W1N9V0iU30LlOzqIB3twRVsxNM83u5wqka29FFUv0w6Fhn+4CqXTvYTO7wDkAkSlypXe3Sn/bVHC90P6u+Aezi/K28bcK9S6i9LEZjT5Vb0O7RMvZPE8qaXS9dL8XibW6ibvggz2U941wt2h2Or+JHdmMl+fC0z8IRG2R1ORcvVHu2zv/ZooS2oDwPv1Vr/NPuC1vrLwPuAj5UiMKcLzluN4Q/KVN8CxNut7j3/ZOneK7aG7OLxGi99lJteLpsTFkXD0kwryubPVaEJahRwaJDX9wE12WeTneoLzqtf5TQxmcFXMvXzL8fw1hFv303izDG7w7GNjD8VV2jR1eByE937Kslwl21xFJqgNgKfUkrlJlUopdzAx7GmoNek7NNreNszmOmUzdE4U7o/Rn/HQTBc+G0ecK1GLl+A+gWXA9Db5rz9fMohGe6iv+OAVd7IQZW4K5m7vong7BVgpm2tPVpogvoEcAtwWCn1qFLqEawW1a1ATY5BAfinKDyjJpAKn8lVUBbni5/ILNBtmY7LV64NmGtLrptv2wZMM21zNOV3rrzRAlkEXkS5B/C2p23bJ6rQHXW3Ym25fjewF9gFfAWYq7XeUrrwnM3a/iA7BiDdfIOJt0v9vVKrm7YQT9M4Uj2niR3cbnc4ZRc98BoAgZky1b6YgnNX4Ao00H/yMP0dB2yJoeBp5lrrTuAxrH2fngE2aK1rvt5/Q+s1gEFUv0yqxqf6Dia7/kkW6JaOYbhyU857a2xWqWma9O3LJiiZXl5MhtubKxpr1wN4odPMQ0qpnwIa+BnWzrk7Mt199aUM0Ok8TeMIzFiMmUoQ2SlrovKZpnmuxNEUSVCllF32ENn1Iul47ZTM7O84SCrShTs0Gl/LdLvDqTq5br4dz2Gmyr8mqtAW1DeAVuAKIADUZf4+CfhaaUKrHKHMl4N0850v2X2SVKTbWqA7aqLd4VQ176gJ1E1dkFkT9aLd4ZRN335rhCE4e5mssSsB34SZeMdNIx3tIbr31bJfv9AE9U7gw1rrl7XWZubPy8BHgd8vXXiVoV5djuGrI370dfo7j9odjmPkxp8mzZUvjzLIVQCooQelaLZ7b/YKmyOpTueNs9uwUWuhCcoFnB7k9TNAqHjhVCaXry63JqqWvhwuRMafyiu0YA2G10/syC4SZ47bHU7JpeNRYu27wXDJ+FMJhRZfDYaL6L5XSUW6y3rtQhPUM8AXlFK+7AtKKT/wecC+SfIOktsOvu13siYqI1vBXArElofLH6B+/hVAbUyW6DuwDdIp/JPn4a6r6aHwkvKERhGcvRzSqbKvibqYdVDXAEeUUo8ppR7DWge1mhpeB5WvbuoCWROVJ93fZ01NlQW6ZZWrEdlW/Wuiornxp+U2R1L9GpZeB0Dv1vKuiSp0HdReYAHwZazyRtuBzwILtNa7Sxde5ZA1UeeLH9ub2UF3hizQLaO66YvwNI0j2XOa2KEddodTMtb0cklQ5RKcuzKzJupQWddEFbxfudb6LPDvJYyl4jW0XsPZ391nrYnqC+MO1O7wXKw9u/+TlJ4pp+yaqK7n7qe3bT2BGa12h1QSic6jJHtO465vwjdhpt3hVL3smqieTb+lt209/gmzynLd4bZ8fwVrt9wL0lqvLlpEFSy7Jqrv4DYiO5+nceWNdodkm9gRq2Et40/l17DESlCRXS+SvvGPcfmDdodUdNF91pTnwKxlGIbsp1oODUuupWfTbwlvf5Yxb/5DDLe35Ne80Jbv4iKFll5H38Ft9Latr9kEZZpp4tkddKV4Z9l5R02gbtpCYod3Et71Ao3Lrrc7pKLLVo8IzpLuvXLxTZiJr2U6/ScPEd2zmfr5l5f8msNt+f7Fkl+9CtWrN3HaF8jtE+UbO8XukMoucaqddDyKu3EsnsaxdodTkxqWXEvs8E56t66vugSV7o8RO7wTMGT/pzIyDIPQkms58+QP6N36dFkS1MVs+f5updQrSqkupdQspdQ9SqlPlDK4SuTy+gktvBKwZrzUoli7dO/ZrX7BFbl9ovo7q2ufqL4DbZipBP5Jc3AHG+0Op6Y0LF5r7RO1bwvJ8NmSX6/QWnwfAP4T+AWQXQu1G/icUupTpQmtcuV2o9xWm2uiZIKE/ax9oqw1UdW2eDy6dzMAwbmrbI6k9rjrmwjOyewTtf2Zkl+v0BbUx4GPaK2/CqQAtNbfBf4Iazt4kcc/WeEdM4lUpCvXV15LzrWgJEHZKfug1LttQ9U8KJlmmuieTYAkKLs0LCnfmqhCE9RsYNMgr78GTCheONXh/DVRtdXNlwx3kTx7AsNbh2+8VJe2U93Uhdbi8d7qWTzef3y/Vb28caxUL7dJcM4K3PVNJE63W+sdS6jQBKWBwUZa343V1ScGCLWuA8NF5PVNpKI9dodTNvFs62nyXAyX2+Zoatv5hT6r40Epkmk91c9ZKQWIbWK4PYQWrwVK/7kqNEF9BrhXKfVtrJl/f6qU+gXwRax6fGIAT8NoArOWQTpZlr5ap8iOP/llgoQjWKWPDCKvv0yqr/L3F5XxJ2fIdvNFdj5HOhEv2XUKLXX0CFbdPT9WmaMbgBhwudb6NyWLrsKdq19VXYPUw5HxJ2fxNI61Kn2nkoR3PGd3OJck2dNJ/4n9GF4/dTMW2x1OTfO1TMM/cTbpeJSofrlk1xmuksRbgUe11mkArfUO4AMli+Tcdf8KuEZr/Y6816YB3wMuB04Cf661/m2pY7lU9XNX4QqE6D95kPiJ/WUrD2KXdCJO/Ph+wKBOtthwjIal19J3YCu9W9fTtOpmu8MZsWzrKTBzCS6P7wJHi1ILLbmO+PF99LY9bW3JUQLDtaB+DRxTSv2zUmpZSa6eJ7Ot/D8C/zTIr+8D2oAxwJ8A9ymlHP9tb3i8hBZl+2qrvxUVP7YX0kl8LdNxyfYHjhGctxpXXT39J/YR7zhodzgjlpu9N0e695wgtOgqDLeXvgPbSIa7SnKN4RLUZODvgTcBryqltimlPqGUmlSSSOBhYCbw7fwXlVLzgFXA57TW/Vrrp4HfAB8qURxFle3mC+94BjOZsDma0rJW90PdtIU2RyLyubx+QousJ9xKrbSfTsTpO7gNgOCclTZHIwDcgRBNl9+Kr2Uqhqc0dfmGK3V0EvhX4F+VUjOAO4C7gK8qpdYDPwIe0Fr3FXKhzGaHowf5lam17gDu0FofU0p9AZiY9/uFwGGtdSTvtd1YY2KO558wE9/4mfR3HCDy+su5KhPVKHZEEpRTNSy9jp7NjxLe9jvGXPe+shT6LKa+A22YyX78E2fjaRhldzgiY/S6Oxi97o6Snb/QSRIHtdZf1VovBZYCL2PN7OtQSn2/wGutAY4P8udo5hpD1WMJAdEBr0WBiinRnL/ZV7UyU8lzFSSmLrA5GjGQb8IsfC3TSPf15qZqV5KIfgmwuitF7bjoOvVa653APcBXgT1YrapC3rdBa20M8udCe1JFgIE73gWB8MXGbpfQ4qvB7aFv/1aSPaftDqck4if2YybieMdMwhNqtjscMYBhGDQsfTMAva9V1oOSmUoSff0VgLIUKBXOcTHFYpuUUu9XSj0MnAA+jVWbb3apgsvYCUxTSuUnqfmZ1yuCO9BAvXoTYNLbtsHucEoiN/40Vbr3nCq06Gpweejb/xrJnk67wylY36EdpGNhvGOn1OTuALVs2ASVl5QeAjqAf8Ta8v1qrfVCrfXdWutDpQxQa62BrcDdSim/Uupa4FbgJ6W8brHld/OZZtrmaIpPJkg4n7u+ifp5qzKFPn9ndzgFi+zeCEC9ktZTrRkyQWVaSh3AN7G62W4HJmmt/0Jr/UqZ4st6F7AAaw3Ud4EPaa23lzmGSxKY0Yq7cSzJro7cl3m1MNMpYkd2ARCQBOVo2QelnteeKnmhz2Iw0ymir1vjT9K9V3uGG/8JAR8F7tdal62YnNb6C4O8dgSo3BWGgOFy07DkWrqeu5/erU8TmF49K+H7Tx0hHY/iaRqHp2mc3eGIYQRmLcMdGk3y7AliR3YSmLbI7pCGFWvXpCLdeJpb8I2fYXc4osyGbEFpra/RWn+vnMmp2mW3P4jsepF0LHKBoyuHdO9VDutBaR1QGZMlct178y+X4rA16KJn8YmR8zaPp25GK2ayv+LrouWTCRKVpWGZNZsvsusFRz8omaaZm15eP/8Km6MRdpAEVWaN2am+W5+yOZLiME1TFuhWGO+oCdRNX2Q9KO183u5whhQ/tpdUz2ncDaPxT5pjdzjCBpKgyiyorLpo8eOVXRctK3HmGKlIN+76ZryjJ174DcIRzq2Jcu6DUkSfm71nGPJVVYvkv3qZubz+c5t9OfjLoVCxQzsAq3qEjBFUjvr5l+PyB4kf30v/yZKuFBkR00wT3v4sAPULpHuvVkmCskH26TW8/RnSyX6bo7k0fQetrcQDsj9PRckvINvjwAel2KEdpHo78TS1UDdV9harVZKgbOCfMBPfhNmkY2GimUHgSmSaafoOWsvR6mYssTkacbGykyXC251XaT+7C3Vo8Vrp3qth8l/eJo3Lzi2YrFT9HYdI9/Xibhwr408VyCogO8MqIPt66XZFvVjpRJxwZnp5qHWtzdEIO0mCskn9oqsxPD5iB7eROHvC7nBGJLs/T2BGq4w/VSDDMGhYfj0Ava89aXM050T3bMKMR/FPnINvzGS7wxE2kgRlE3ddfW7wt1InS5wbf2q1ORIxUqHMg1LfgTbHPCiFt1l1AqX1JCRB2Sg7BtDbtgEznbI5motjphLEDmfq70mCqljuQChX484J+5WlIt1E978GhovQwqvsDkfYTBKUjeqmLsQ7ehKp8Bmie1+1O5yLEj+2FzMRwzt2Cp6GwTZKFpUi1823db3tD0rhnc9DOmXVDKxvsjUWYT9JUDYyDONcK2rLEzZHc3H6DpwbfxKVzUkPStnZew2t19gah3AGSVA2a1hyLbg8RPdtqaxN5GT8qWqc96Bk42SJ+PH9xI/tweUPEpx3mW1xCOeQBGUzaxO5y8BMO2IMoBDp/hixo3vAcFFXRduG1LJQ6zpwuYnufdW2B6WezY9YsSy9DpfXb0sMwlkkQTnAuTGApypit93YkV2QTuKfMAt3Xb3d4Ygi8ISaCc5dlXlQKv+s0lRfb67Cf9PKG8t+feFMkqAcIDBzCZ6mFpLdp+g70GZ3OBeU696bKd171aRplbUnaM+rj2OmkmW9du/WpzGT/QRmLcM7elJZry2cSxKUAxiGq6ImS0T3bQEgMHOpzZGIYqqbvhjv2Cmkwmdz+zCVg2mm6dn8GACNK28q23WF80mCcoiGJdeC4SLy+iskw112hzOkRPdJEqeOYPgCUsSzyhiGQePKTCtq0yNlu27fvi0kuzrwNI0jOGdF2a4rnE8SlEN4GscQnLsS0inCbc6dLNGXmYYcnLUUw+21ORpRbA2t12D4AsSO7CrbfmXdmWTYuOJGDJe7LNcUlUESlIM0Lr8BgJ4tTzp2skR2nUxgtjzpViOXP2C15ilPK6q/8yh9+17DcHtz3dxCZEmCcpDArGV4GseS7OrILYR1knQinisQG5QEVbUaM7PowtufIdUXLum1zj77M8Ak1HoN7mBjSa8lKo8kKAcxXG4almWmnG953OZo3ih2eCdmsh/fhFl4GkbZHY4oEd/YKQRmLsVM9pd0ynn/ycNEdjwPbg+jrrq9ZNcRlUsSlMM0LL0ub7LEWbvDOU9072YAGciuAU2rfw+Arhd/RTreV5JrnH32p4BJ47Lr8TSNK8k1RGWTBOUw+ZMlereutzucHNM0c+NPwTkrbY5GlFpg9gr8UxTpaA/dLz9Y9PPHTxwgsnsjhttL85p3Fv38ojpIgnKgxuVvAaD3tSccM1ki0XmUZFcHrmAj/omz7Q5HlJhhGIy+9n0AdG38NalId1HPf/aZ+wBoWHkjnsYxRT23qB6SoBwoMGspnqZxJLtO0rd/q93hABDdl51evkymAteIwLSFBGavwOyPcfaFXxTtvLGjrxPdswnD66f5ituKdl5RfSRBOZDhctOQnXL+6mM2R2OR7r3aNPraOwGDns2Pkug+ecnnM5MJTj38TQAaV92MJ9R8yecU1UsSlEM1LH2zVV16z2bbt+FIRXuJHdoBLjeBWVLeqJb4x88gtPhqSCU5u+F/Lvl8Z5/9GYlTh/GMmsCoq36/CBGKaiYJyqE8oWbq1ZvATNNjc32+iH4JzDSBGUtwBxpsjUWU36i17wG3h/D2Z3IVx0cidvR1ul78FWDQcsuf4/LVFS9IUZUkQTlYdsFk72tPlr26dL7IrhcACC1cY1sMwj7eURMYe8MfAXDqt9+kv/PoRZ8jnYhz6jf/BmaapsvfLnUcRUEkQTlY3bRFeMdMJhU+S3TPJltiSEW6reoRLg/BeattiUHYr2HFjdQvvBKzP0bHA/eQTsQLfq9ppul87LskzhzDO3YKo655bwkjFdVEEpSDGYZB4wprynnPq/ZUlsh1781sxR0I2RKDsJ9hGIx760fwjp5E4tRhTj/6HUzTvOD7zHSKUw/9B71bn8Zwe62uPY+vDBGLaiAJyuFCreswPD76DmwlceZY2a+f695bIN17tc7lDzD+XZ/A8PgIt22g42dfJRXtHfJ4M5Xg5K/+mXDbBgyvnwnv+Qz+SXPKGLGodJKgHM4dCFG/8Cqg/K2oVKSbvkM7pHtP5PhapjP+XZ/EVVdPdO9m2r/7cWJHdp93jGmmie7bwvGffJnIrhcx/EEm3vE5AjOX2BS1qFQeuwPIUkp9DPgYMAbQwMe11s9mfjcN+B5wOXAS+HOt9W/tirXcmlbdRLjtaXq3Ps2oa+7A5fWX5bqR3RvBTBOcvVy690ROcM4KJv/xPZz85T8TP/o6x370WbyjJ+AdOwVP4ziiezeR7LLWTLkCISbe8TmpPiJGxBEtKKXUO4G/Ad4GjAK+CTyklMpWkLwPaMNKXn8C3KeUmmVHrHbwT5yNf9Jc0rHIJU3zvVjhTPdevXTviQG8TS1MuuvLNF3xDnC5SZw5TvT1V+jZ9FuSXSfxNI1j1Lo/YOqH/1WSkxgxp7SgJgJ/r7Xemfn5+0qpfwJalVLtwCrgBq11P/C0Uuo3wIeAv7Mn3PJrXHkTp47toWfzozQsvQ7DMEp6vWTvWWKHd4LbQ3DeZSW9lqhMhtvDmOvuYvTa95I4c4z+0+0kzhzHP34mgdlSEktcurIlKKWUDxg9yK9MrfV/DDh2LRACdgBXAIe11pG8Q3YDNTUoUr9wDZ1P/oD+E/uJH9tD3eR5Jb1e79anrO69OZfhrqsv6bVEZTM8Xnwt0/G1TLc7FFFlytnFtwY4Psif81b9KaUWAz8FPqu17sBKVNEB54oCwVIH7CQujy+3JXapt+I20yl6M9UrstPchRCi3MrWgtJabwCG7ZdSSr0N+DHwD1rrr2dejgCBAYcGgdLuRe1AjSveQveLvya86wXGXP8B3PVNJblOdN8Wkj2n8TSPl5lXQgjbOGKSBORm8f0P8Mda66/l/WonME0plZ+k5mderyne5vHWbrapJD2vlW4r7t7MdPbGFW/BMBzzERFC1BhHfPsopd4N/D1wvdb6gfzfaa01sBW4WynlV0pdC9wK/KT8kdqvcdXNAPRsfgQzlSj6+RPdJ62tNdweGpZcW/TzCyFEoZwyi+9TgB94SimV//p7tdYPAe8CvoO1Buo08CGt9fayR+kAgVnL8I6dQuJ0O+FdL9KweG1Rz9+75UnAJDT/ipJ1IQohRCEckaC01isu8PsjwM1lCsfRDMOg6U23cPrhb9K98TeEFl1dtCnnZipBb6brsEEmRwghbOaILj5xcUKL1+Kub6K/4wCxwzuKdt7I7pdIRbrwjp1C3dQFRTuvEEKMhCSoCuTy+GhceRMA3Rt/U5RzmukUZ5/9KQBNq95a8oXAQghxIZKgKlTjihsx3F6iezePaAO5gXq3Pk2i8xieURNy662EEMJOkqAqlLu+iVDrNQB0v/TQJZ0rnYhz9pmfATD6mjsw3I4YmhRC1DhJUBWs6U23AAa9W5+m/9ThEZ+nZ9MjpMJn8I2fSb1s6y6EcAhJUBXMN3YKDStugHSSU7/9NqaZvuhzpPrCdL3wCwBGX/c+WZgrhHAM+TaqcKOvfR/u+mbi7btzU8Qvxtnn7icdi1A3o5XAzKUliFAIIUZGElSFc9fVM+YtHwTgzNM/Jhk+W/B7w7tepOflh8BwMeba98nMPSGEo0iCqgL1C9YQmL2cdCxC5xPfL+g98Y6DnHrw3wAY/ea78E+aU8oQhRDiokmCqgKGYTD2pj/B8PiI7Hyezqd+OOx4VCraS8fPv4aZiBNqvYam1beUMVohhCiMJKgq4W0ez7i3/Rm43HRv/A0nf/kN0sn+NxzXf+oIJ3729yS7TuKfOJuxN39YuvaEEI4kC16qSGjR1biCjXQ8cA+RXS+S7D1D44q34B09CXewka4Xf03va0+CmcYdGsX42/8Wl9dvd9hCCDEoSVBVJjhzKZP/8Cscv+8rxNs1p9r1+QcYLhpX3Miote+RauVCCEeTBFWFfC3TmfxHX6dny+MkTreT6DxGsvskddMWMvq6u/CNnWJ3iEIIcUGSoKqUp2EUo9e+x+4whBBixGSShBBCCEeSBCWEEMKRJEEJIYRwJElQQgghHEkSlBBCCEeSBCWEEMKRJEEJIYRwpFpYB+UGOHHihN1xCCGEGCDvu9k98He1kKAmAtx55512xyGEEGJoE4F9+S/UQoJ6BbgaOA6kbI5FCCHE+dxYyemVgb8wTNMsfzhCCCHEBcgkCSGEEI4kCUoIIYQjSYISQgjhSJKghBBCOJIkKCGEEI4kCUoIIYQjSYISQgjhSJKghBBCOFItVJIoiFJqKfAtYAmwH/ig1voNK5uVUtOA7wGXAyeBP9da/7ZMMd4A/AMwN3Ptf9Raf3uQ464DngD68l7+mtb6y2WI8YPAt4F43ssf1Vr/cMBxttxHpdSdmfjyBYCntNZvGXBs2e+jUmo18JDWuiXzsw/4d+B2rEoo39Baf3WI9xrAl4H/BfiA7wOf1FonSxxjC/AvwJsBA3gE+JjW+uwQ7z8MjAGyVQKOaq1ViWP0A71Af95hLwz8b5451q77GB5wiAfwA5O11scGeX/J7uNQ3zXl/jxKgiL3JfBr4F5gLfAu4HGl1HStdc+Aw+8DXgR+D7gK+JVSapnWen+JY5wKPAC8PxPrSuAxpdRBrfVjAw5fAdyvtX5vKWMawgrgn7TWn7rAcbbcR631fwP/nf1ZKbUceBz45CCHl+0+Zv5n/hBwz4BffRFQwGygCXhUKXVUa/2jQU7zv4B3YsUdB34JfAb4Uolj/C7QDcwEvMCPgf8A/mCQc4wFJgONWutIMeIqMMZW4IzWekIBp7HlPmqtQ3nHeID1wIYhklPJ7uNw3zXAOsr4eZQuPss6wKu1vldrndBa3wfsAN6Tf5BSah6wCvic1rpfa/008BusD1upzQB+orX+pdY6nWndbQCuHOTYlcBrZYhpMBe8ts33MT8OL1ay+oLWeusgh5TzPn4R+AjwlQGvvx+4W2t9Vmt9EOtL7cNDnOP9wL1a63at9SngC8McW5QYlVIuIA18UWsd0Vp3Af+F9dAxmJXAnlIkp6FizLtuof8ty34fB/G3WMn+80P8vpT3cQZDf9eU9fMoLSjLQmDXgNd2Yz11DTzu8IAPxW5gdQljA0Br/SzwbPZnpdRorCK4Px7k8BXAOKXUR7C6XH4KfFZrHR/k2KJRSrmxukjvUkp9A4hiPV1/TWudX/TRtvs4wEexuu/+c4jfl/M+fktr/Tml1LrsC0qpZqwimjvzjhvsc5m1cJBjJymlRmutz5QiRq11GnjHgOPeAWwZ4hwrAJdS6mWsFterwF9qrQf+/1e0GPOu26KUagPGA89krnt0kHOU/T7mU0pNwmppXJm5v4Mp2X28wHdNWT+P0oKyhLC+TPNFgeAIjysppVQTVovjJawmeP7vPEA7VnN6AXAdcD1WX3CpjQM2AT/E+p/mdqwnxY8MOM72+5jp1v0kVuvpDRWTy30fB+vGwbpPcP69Gu4+Dbyv2b8X5b4OEeN5lFKfwEpQfzvEISngZayun+lYiewRpVSpY4wAz2ONkymsB5NfDnGs3ffxr4BHtdbDtfhKeh+zBnzXbM68XLbPo7SgLBGsgfJ8QWDgoGWhx5VMpnvs11hPJncOfMLKDEC+Oe+lvUqpu4GvAX9Tyti01ieAa/Jeek0p9W9YY3r5rRTb7yNwE1bX1MOD/dLO+5gn28LMv1fD3aeB9zX7RVDy+5rpLv034BbgOq317sGO01p/fcD7Pg38GVaX1bODvacYtNZ/PeC6fw2cUkpN1VofGXC4nffRjdU1NuwGduW4jwO/azh3T8r2eZQWlGUn1lNVvvmc3zzNHjdNKRW4wHEloZRai/Uk8yvgdq11bJBjJiulWc9tyQAABWBJREFU7sm0ELJ8wBuOLUF8i5RSXxzw8mDXtvU+ZtwK/GyoLhQ772NWZhbcCc7/bA53nwZ+jucDxzPjQiWjlGrAmu14GbB6uCd/pdRfKqXyx6fcWA/KJb2vSqkvKaUW5L2U/e862HVtuY8ZazL/fGq4g0p9Hwf7rrHj8ygtKMt6wFBK/RXWFMp3YY2lnNcFoLXWSqmtwN2ZJ5Y1WF90V5Q6QKXUbOAh4O+01v82zKGdWE87UaXUl7C62j4L/N9Sxwh0AR9XSrVjTSFfDvwF8L/zD7LzPua5HPg/w/zezvuY78fA5zNjJyHgE1hTuoc69hNKqaewnl6/wOBjlMV2H9bD7tVa64FdtwPNwBqjfBvW5+VrwB6sMZRSWgKsUkplZxb+C/BwZvB+ILvuI1ify43DjD1lzaBE9/EC3zVl/TxKCwrQWvcDN2MlpjPA3wHv0FqfUkrdOWB9wruwxiROYk0A+JDWensZwvwo0AB8VSkVzvvztfwYM62qm7Gmy3diDQbfD3yj1AFmBpzfjjVTpwdrquqXtdY/d9B9zJoBnDcW4JT7OMDngO1Ys0pfwbqn38qLOZxZ20Xm9fuBF7C+rHZm3l8ySqklwFuxJriczPtctg8R46eAjVhjJieBWcAtWutS73b9IeAssBc4iLUe6q4hYiz7fcwzgwGfy6wy3schv2so8+dRdtQVQgjhSNKCEkII4UiSoIQQQjiSJCghhBCOJAlKCCGEI0mCEkII4UiSoIQQQjiSLNQVokiUUj/AKlMzlC9iVYVeDzRorctS2ilTPud54A+11q8Pc5wLa23NXVprXY7YhBiOtKCEKJ6PYVV7noi1hQtYC1izr92DtWhxIufq7JXDXwBbh0tOkKtM/iXyFl4KYSdZqCtECSilFgPbgJmZfXPsiqMOOIxVwLWgSh1KqX1YlT02lDI2IS5EuviEKKPMHkC5Lj6llAncAXwaq7DmJuB9WFuB3IVVMurTWusfZ97fAPwT1lYmJvA01vbqQ23h8F6gKz85KaX+D9Zup+Ow9kH7jNb6kbz3/BKrNbihCP/KQoyYdPEJYb9/AP4Sq1DoNKyCnz1Y1cF/AXxbKZXdG+o7WInsRqytTUys7biHetj8PeDR7A9Kqdsy13ofVnXph4H7lVKNee95FLh+mHMKURaSoISw339orddntql4CGu/nM9kJip8A2tPnZlKqVlYLaI/0Fq/kmkV3YVVYPSmIc69CquwZ9YMIA4cynQ9fglr07tE3jE7sSpVzy/Kv50QIyRPSELYb2/e36PAwbxdfrP7+/ixdk4F0Eqdt31ZEKtV9dAg5x4PnM77+f9hzTTcr5TajLVb6ve11n15x3Rm/tlykf8eQhSVtKCEsF9iwM9D7QXkyRy7HFiW92ce8P0h3pMGjOwPmf2PVmK1uF4APgC0ZSZ1ZGW/F0q9BYYQw5IEJUTl2AV4gXqt9V6t9V7gOPCPWElqMCewJkMAoJR6J/BhrfXjWuuPYbW8erH2dMoal/deIWwjXXxCVIjMTsS/AX6klPoocAq4G2tyxe4h3rYZWJr3sxv4R6VUB9aMwcuBCZm/Zy3l3OZ+QthGWlBCVJb3YyWTX2HtaNoE3KC17hri+IexZvsBoLW+H/g8VqvrdeArwP/WWj+d9561wKNl2OVWiGHJQl0hqphSKoi1xflNWutXCzjeBRzCmin4bInDE2JY0oISoopprf9/u3ZswyAMAFH0tqFjFdZAyggMkQ0zBjUFqdxAhU7ovd6Su6+zvOdcS5+bR5YkP3GigUDB+32TzNPwN330X09bkvWRW8EFT3wAVLKgAKgkUABUEigAKgkUAJUECoBKBwCM+k7ps1wSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_velocity(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although we compute acceleration inside the slope function, we don't get acceleration as a result from `run_ode_solver`.\n",
    "\n",
    "We can approximate it by computing the numerical derivative of `ys`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXybV5no8Z8kW9733bEd20l8sqclbRra0o1hH6C03DvQUgqUYRmWmQHmA1zmDlOWYe6wDAxlh8vOhYFSKC2Uli7QBtqmabPHJ3G87/u+abt/vHplxXUSW5b0vpKe7+fjj+3XlvTEVfXonPOc5zgCgQBCCCGE3TitDkAIIYRYiSQoIYQQtiQJSgghhC1JghJCCGFLaVYHYDWlVAZwOdAH+CwORwghUo0LqAIOaq0Xwn+Q8gkKIzk9bnUQQgiR4l4EPBF+QRKUMXLixz/+MZWVlVbHIoQQKaW/v59bb70Vgq/F4SRBBaf1KisrqampsToWIYRIVc9bYpEiCSGEELYkCUoIIYQtSYISQghhS5KghBBC2JIkKCGEELYkCUoIIYQtSYIStuD3B+genKKzf9LqUIQQNiH7oISl/vB0J48f6eF0xxjTcx4ArrlkA+943S4KcjMsjk4IYSVJUMISgUCAH/7uFD9/+EzoWnF+BtNzXv50uIfDZ4Z45+t28aJLNuBwOCyMVAhhFUlQIu4CgQDfvvc49/6pFafTwdtfs5P9O6soLcykf2SWu35+mKMtw3z2R4eYmffyihfWWx2yEMICsgYl4ioQCPC1u49y759aSXM5+MibL+PVL2qkrCgLh8NBVWkOn3rXldzxmp0AfO++E4xOzlsctRDCCpKgRFw9eqib3/2lHXeak39+2xW8cFf1837H4XDw2msauXx7BbPzXr796+PxD1QIYTlJUCJuZuY8fPe+EwC8++Y97N1acd7fdTgcvOt1u8lwu3j8cA/PNg/GK0whhE1IghJx85PfNzM+tcC2+mJuuKz2or9fXpzNG1+iAPjaL4+w4JHzJIVIJZKgRFy09U5w3xOtOB3w7pt343SurjLvtdduor4qn/6RWX75yJmL30AIkTQkQYmYMwsj/AF41dWNNFQXrPq2aS4nf3ujUTDx27+04/X5YxSlEMJuJEGJmDtwtJdT7aMU5mVwy8u2rvn2uzaVUluRy/jUAgdP9scgQiGEHSXcPiil1D7gPq11efD7DGAKWAz7tT9rrV9qRXzi+X79x7MAvOElitys9DXf3uFw8NIr6vnOvcd54MmOFSv/hBDJJ2ESlFLKAdwBfG7Zj3YBo1rryvhHJS6mpWuc5o4xcjLTVlUYcT7X763h+/ef5Dk9yODoLOXF2VGMUghhR4k0xXcn8G7gU8uu7wUOxz8csRq/eaIVgL/at5GsjMjfDxXkZnDl7ioCAXjo6c5ohSeEsLFESlBf11rvBZ5Zdv0FQLlS6qhSakAp9XOl1AYL4hPLTEwv8PjhHhwOeNVVDeu+v5ft3wjAH57uwCfFEkIkvYRJUFrr3vP8aAY4ALwYUMAccE+84hLn9+BTHXi8fvZuraCqNGfd97drUylVpTkMT8xzSMvGXSGSXcIkqPPRWn9Aa/0+rfWQ1noc+ABwuVIq8gUPsW4+n5/f/rkdgFdf3RiV+3Q4HLzsCmMU9eCTHVG5TyGEfSV8glJKfUIptS3skjv4WTqMWujJE/0Mj8+xoSyHS5rKona/N1xWi8MBh5oHmVvwRu1+hRD2k/AJCtgNfF4pVaiUKgS+BNyvtR6yOK6U9odgIcMrr2pYddeI1SjKz6Sprgivz8+RM/KfWIhklgwJ6g5gDGgB2jH2Q91mZUCpbnrOw+HTgzgdcM0lNVG//8u3G01mnz4hm3aFSGYJsw/KpLV+DCgM+34EuNWygMTzPH2iD68vwO7NpRTmRf/Y9n3bK/nR75p55tQAfn8gqiM0IYR9JMMIStjME0eMgsur9sSm40N9VT6lhVmMTS1wtmc8Jo8hhLCeJCgRVTNzHp7TQzgc8MJdVTF5DIfDETbNNxCTxxBCWE8SlIiqp0704/X52dFYQlFeZsweZ992o7PVwVOyDiVEspIEJaLqQHB67+rdsW3ountzKRluF2e7JxiZmIvpY4nE5fP56eyf5PHnejjZNkIgELA6JLEGCVckIexrZs7Ds3oQhwOujHGCcqe7uGRLGU+d6OeZUwO8bH99TB9PJJaugSm+dvdRmjtG8XiX2mLt3lzKba/cxtaNxRZGJ1ZLRlAiap4+aUzvbW8ooSg/dtN7JlmHEit56ngfH/zSnzh2dhiP1095cTaXbasgJyudoy3D/NN/Pc6/fe9p5hdlo7fdyQhKRE1oei9G1XvLXbbNSFCHzwyx6PHhTnfF5XGFPfn9AX72h9P85PfNgFFF+p7X7yEv22guMz3n4Z7HWrj3T2f5y7E+vnb3Uf7hDZficMg2BbuSEZSICo/Xx+FgZ4f9O2NTvbdcSUEW9VX5LHp8nO4ci8tjCvv62UOan/y+GYcDbn/Vdj5822Wh5ASQm5XOba/Yxmfffw3udBePPNPFA9LT0dYkQYmoaG4fY2HRR11lHqWFWXF73J2NJQAcbx2J22MK+znZNsJPH9I4HPC/3rKP19+w5bwjo/qqfN73P/YA8M17jsmbGxuTBCWi4rnTxvEXlzaVx/Vxd24qBeD42eG4Pq6wj+k5D5//8SH8Abj5+i2rGsFft7eWV13VgNfn5zPfP8jkzGIcIhVrJQlKRMVzwfOZLlXR61y+GjuCI6hT7WPnVGuJ1BAIBPjqL44wODbHltpCbn351lXf9o7X7ERtLGJ4fI7//sPpGEYpIiUJSqzbxPQCZ3smSHM5QwkjXgrzMqityGXR46OlS9oepZpHD3Xz+OEeMt0uPvSmvaS5Vv+Slp7m5N037Qbgd39pZ2xKTuixG0lQYt2OnBkiEIAdjcVkuuNfGLqzMTjN1yrTfKlk0ePj+/efBOAdN+6iujR3zfexqaaQK3ZUsujx8ctHW6IdolgnSVBi3Z4LHr0V7/Un085NUiiRih54sp3RyXkaqwt48eV1Ed/PG16qAPjtn9sZn1qIVngiCiRBiXUJBAIcNgsklDUJKrQO1TaCzyfrUKlgftHLzx8+A8AtL1PrOnJlc00h+7YHR1GPySjKTiRBiXXpHpxmeGKewtwM6qvyLYmhpCCLqtIc5hZ8nO2ZsCQGEV+/C452NtcUsG9H5brv742hUVSbjKJsRBKUWBezeu+SpjJLDw4M7Yc6K9N8yW5uwcvdjxqjp1tfvi0qnSA21xZy+fYKFhZ93Pv42XXfn4gOSVBiXZ47HVx/inN5+XKh/VBSKJH07j/QxsT0IqquiL1bozetfPP1WwB4+GAXPr90PbcDSVAiYh6vn2PBDbKXWFQgYTILJU62jsiLSxLzeP38+o/GCOeWl2+Nah+97Q3FVJXkMDo5z5HgGy9hLUlQImJne8ZZWPSxoSyX4jh0L7+Q8qJsyouzmZn30t4r61DJ6ukT/YxPL7CxMo9Lm6I7anc4HNxweS0ADz/TGdX7FpGRBCUidrJ1FCDum3PPZ3uDccaPlt5qSeuBJ9sBeOn+jTHpQn79XiNBPXmsj5k5T9TvX6yNJCgRsZNtRkGCmRispuqKANAdkqCSUf/IDIdPD+FOc4YSSbRVFGeze3Mpi14/TxzpicljiNWTBCUiEggEONlmjKC2N9hjBNUUTFDSnTo5PfiUcTTGVXuqzzlGI9puuCw4zXewK2aPIVZHEpSISPfgNFOzixTnZ1BZkm11OAA0VOeT5nLSPTjNtEzPJBWvz8/DB411oZftr4/pY125u5pMt4tT7aP0Dk3H9LHEhUWcoJRSTqXUdqXUdUqpFymltiil5GjKFGGOnrY1lNjmRNL0NBebNhQAcEZGUUnl4MkBRicXqCnPjfmUclZGGlfuNk6FfvgZGUVZac0JSil1jVLqbmAMOA48AvwRaAaGlFI/UkpdGd0whd3Ybf3J1LRRpvmSkTm997L99XF5Q/TiYDXfn57rJhCQbQtWWXWCCo6QHgH+L9AK3ARsADKBbKAeeAvQC/xUKfWoUqop2gELezhls/Unk7kOJZV8yWN0cp5DzQOkuZxcv7cmLo+5o7GU/Bw3/SOzdA5MxeUxxfOt5WyEHwGf0Frff56fdwU/7lNKfRi4MXibfesLUdjN6OQ8fSMzZGW4aLCo/975qLBCiUAgYJvpRxG5vxztJRCAvVvLKcjNiMtjupwOLt9ewcMHu3jqeD8bK+31PE8Va5ni23+B5HQOrXVAa30PcEVkYZ2fUmqfUmow7Hu3UuqbSqlRpdSQUuqj0X5McS5zek9tLMa1hgPi4qGyJJv8HDcT04sMjs1ZHY6IggNH+wC4ek91XB/3ih3G0fFPn+iP6+OKJat6dVFKKeAGpVTesut/faHbaa2jNnmrlHIopd4OPAiE15jeCShgE3A5cLtS6s3RelzxfHYrLw/ncDiWys1lP1TCG5ua50TrMGkuJ5dvX3/X8rW4tKkMd5oT3TnG6KSctmuFiyYopdR7gF8D7wdOKKVeF/bjT8QqsBXcCbwb+NSy67cDn9Zaj2mt24HPAe+MY1wpxxxB7Wi0V4GEqam2EJB1qGTw5LE+/AGjGXFOVnpcHzszI409wXZKMoqyxmpGUO8E9mqtXwtcC/yzUuoDwZ/Fc4L/61rrvcAz5gWlVCFQBZwM+71mYFcc40ops/Me2nomcDkdNNUWWR3OiqSSL3kcONoLwFW74zu9ZzKn+Z6SBGWJ1SSoNK31DIDWug24DniZUuoLxDFBaa17V7icG/w8G3ZtFqOqUMTAmc5x/AFo3FBAZsZaamzix5ziO9s9jldO2E1YE9MLHDs7gsvp4IooHEoYiX07KnA44MiZIeYWvJbEkMpWk6D6lVKXmN9oraeAVwGlWD9SmQl+zgq7lg3I9u8YOd1ljErMajk7yst2U12aw6LXT3vvpNXhiAg9ebwfvz/AnqYycmPY2uhCivIyaaorwuP1hw7nFPGzmgT1ZuCc8a3W2qu1fjNwTUyiWiWt9RhGbCrs8lbOnfITUXSmaxyALXWFFkdyYeY0n6xDJa4DwWatVk3vmczRm0zzxd9FE5TWultrveJ/Ga31geiHtGY/BD6ulCpVStUDHwpeEzFgthDaYtP1J9OWYKHE2e5xiyMRkZicWeRIyzBOC6f3TPt3GutQB08OyGGYcRbRJpZgMrDLBph/wWi5dAI4CNwNfN3SiJLU2OQ8wxPzZGWksaEs9+I3sNCmDWaCksMLE9Gh5gH8/gC7NpXEbXPu+dSU51JRnM3U7CKtPfKGJ54iXeX+BfAb4PMQ2g+1H7hHa30oSrGtSGv9GFAY9v088J7gh4ghc3pvc00hTqe9OzQ0bijA4YCO/kk8Xh/paS6rQxJrcOiUsd4T771PK3E4HFzSVMbvn+zg8Okh288eJJNIR0G7gD8AKKVqMUYtrwH+qJS6LjqhCbsxCySabL7+BIRGeT5/gI4+6aWWSHz+AM8GCxJeoMotjsZwSXA/1OHTQxZHkloiTVCZwGjw69cDT2mtdwP/jDHlJpLQmc5ggUSCvIMMTfPJtExCOds9ztTsIuXF2dSU22MqeffmMhwOo4vK/KKUm8dLpAmqFdge/Pom4KfBr+8F9qw3KGE/gUCAM11mgYT9R1AAm2qMs6FkHSqxHGo2Rk97Vbltmv3m57jZtKEAr88favUlYi/SBPVN4GtKqW9jNIT9TfB6Fuf2yRNJYmB0lqlZDwW5bsqKsi5+AxsIJSgZQSWUQ80DgNG93E72bDGm+Y7INF/cRJSgtNZfBr4AVALv1Vqbx05eAXREKTZhI+HTe3Z5V3sxjcEpvrbeSekokSAmZxY50zlGmsvBrs2lVodzDlmHir+Ie9Vore8C7lp2uYSl6T6RREIFEgkyvQeQm5VOZUk2/SOzdA1M0VBdYHVI4iIOnx7EH4CdDSVkZ8a3OezFbG8owZ3mpLV3gonpBcvL31PBqhOUUuoMxlEXDwGPaK2f10NGa/3ZKMYmbGSpg0RiFEiYNtUU0j8yy9nuCUlQCSC0/mSz6T0Ad7qL7Q0lHD4zxJEzQ1xzaXxO901la5niexVGC6G3Au1KqQNKqY8rpa600aZdEQM+fyDUkSFRCiRMmzbIOlSi8IeVl+/dWmFxNCuTab74WvUISmt9GjgNfEUp5QJeCLwEY7PuFqXU4wRHWFrrllgEK6zRPTDF/KKP8qKshJvW2FQjHSUSRVvvBONTC5QUZFJXmXfxG1hgT1MZ3A+HzwwRCAQSZj02UUW0BqW19gFPBD8+rpQqAF6MkbA+AGyJWoTCcqHpvQTZ/xTOHEG19U7g8wdw2bwDRioL35xr1xf+xuoC8rLdDI3N0Tc8Q7XNW34luoiLJJRS2UADYL6lbge+FfwQScScHjPLthNJQW4GpYVZDI/P0Ts0TW2FPd+ZCzjaMgzApU32W38yOZ0Odm8u5cDRXo6dHZEEFWORNot9EzAIHMU44fYZjEat5meRRNqCZyqZnRkSzdI6lEzz2ZXH6+dUu7EBduemEoujubAdjUZ8J1qHLY4k+UVa3PAZ4DtAI8aR61VAddhnkST8/gCtwRf2hg35FkcTmc1y9IbtnekaY2HRR21FLkX5mVaHc0FmAj3ROmJxJMkv0im+fOAurbVsyk1yA6OzzC14Kc7PoCjP3i8c59MYHEG1ygjKto4Fp/d2bbLX5tyV1FXmk5OVzuDYHIOjs5QXZ1sdUtKKdAT1Q+AtUYxD2FRrb3D0lMB7iBqrlwolAgE5cM6Ojp0NJiibdY9YicvpYEdDcBTVJqOoWIp0BPVZ4Fml1K0YxRHn9JHRWt+wzriETbQFRx3mKCQRlRRkkpedztSsh5GJeUoLE6OXYKrweH2cCjZg3dlo/wQFxjrU0yf7OX52hOv31lodTtKKNEH9EJgG7gdmoxeOsJtkGEE5HA4aqgs42jJMa++EJCibOd05zqLXT11lHoV5ibHPbmkdSgolYinSBHU5cIXW+mg0gxH2Y67bbErgERRAfXU+R1uGaeudYJ8NTmkVS8zy8t0JsP5katxQQKbbRc/QDGOT87Yv7EhUka5BacKOXRfJaWJ6gZGJeTLdLipLcqwOZ12W1qGe10JSWCxUIJEA60+mNJeTrfXFgKxDxVKkI6jPAN9TSt0FnAU84T/UWv92vYEJ67WFTe85E7wDgzlF2SaVfLay6PHR3GGsP5n7ixLFzk0lHD49xImzI1y9Z4PV4SSlSBPU/wt+/twKPwsArgjvV9hIaP9TdWLufwpXW5GLy+mgb2SG+QUvmRkRN1ERUaQ7xvB4/dRX5Sdcn0ezoOO47IeKmbUct1GotR4H0FqvampQKVWktR6LNDhhrdYeYzqsMUE7SIRLT3NRW5FHe98k7f2TbN1YbHVIgqXy8t0JNL1n2lJbSHqak/a+SaZmF8nLlsPEo20ta1B/VEp9ONgY9oKUUqVKqX8G/hR5aMJqZgVfY4J2kFjOHAnKNJ99HD9rjD7s3t5oJe50F03B89Gkq0RsrGWe4yrgU0C3UuoA8ABwAhgGHEAZsAe4FngR8P3gbUQCWvD46Bmcwul0UFeZLAmqgEcPdUuhhE14fX50pzHBsr0h8RIUGOtmJ1pHONU2yv6dVVaHk3TWch7UNPAPSql/B94J3AJcytJ6kwd4DmNv1N9qrXujHKuIo46+SfwBqKvIJSM9OZYUwztKCOu19kyw6PGxoSwn4dafTNuClXxmo1sRXWteKdZa9wN3AncGT9ItAfxaaxnjJpHWJOggsVx9cIqvvW8Svz+Q8JWJic58Ud9Wn5ijJ4CtG40pvpbucTxeH+lpyfFmzi7WVcqktfYDcvZxEgqtPyVwB4nlCnIzKM7PZHRynv4ROWzOamZ7I3M/USLKzXZTV5lHZ/8UZ7snEvrfYkeRbtQVSa49uE5TX5Uc60+mUKGErENZKhAIcKrdmHTZ3pDYL+rmNN/JNpnmi7akSVBKqbcppTxKqemwj9utjisRBQIBOvuTM0E1bpB1KDsYGJ1ldHKBvOx0NiT4SNZMUOaGYxE9ybRb8QXA57XWH7E6kEQ3ND7HzLyX/Bx3wjTvXK2GquDZUJKgLGWuP22tL074tcBQoUTbKIFAAIcjsf89drKuBKWUSg/exzn/RbTWVnQ43wt8yYLHTTqd/VOAMXpKtv/Z6mWKzxbM9adtSbBmU1WaQ0Gum/HpBfpHZqkqTey+lXYS0RSfUmq/UuoIMI9x7MbUso+4Ukq5gN3AbUqpXqVUi1LqI0qp5Hp1jZP2PuPFe2OSTe8BVJfl4k5zMjw+x/TsotXhpCxzBJWo+5/CORyOUGcSc11NREeka1BfBCaAG4EbVviItzLgGYzNwQ3A64F3Bz/EGnWYCSpJNuiGczkd1FUtlZuL+Jue89DRP0may8Hm2sRvowXh+6Gks1s0RTrFtwvYr7U+Fs1gIhXcm3Vt2KXDSqkvAzcDX7UmqsS1NILKsziS2Gioyqela5z2vkl2JtAZRMlCd4wSCMCmmsKk2QS+rcFch5IRVDRFOoI6BVRHM5D1UErtUErdueyyG2MKUqyB1+ene9CYpa2rSM4EVS8jKEsl0/qTaXNNIWkuJ50DU0zPeS5+A7EqkY6gvgx8KzhKOQOcM5lvwXlQ48AHlVLdwHcwWjC9H3hvnONIeL1D03h9ASqKs8nOTLc6nJgIdZSQQglLLK0/JU+Ccqe72FxTQHPHGLpjlL1bK6wOKSlEmqC+G/z8f1b4WdzPg9Ja9yilXgP8B/CfGA1sP6m1/kU840gGHX1LFXzJylxba++Xlkfx5vMHOB1sEJtsXRe21hfT3DHGqTZJUNESUYJa7XlQ8aS1fgS4zOo4El17f/JW8JnOaXk0OkN1aWJvFE0knf2TzC/6KC/Opigv0+pwompbfTG/+uNZ2bAbRevdB/ViYAfGWtYp4GGttTcagQlrmBV89UlYwReuvjqf0cl52nsnJUHFkTl6UsFzlJKJCjaOPd05js8fwCUj83WLdB9UpVLqSeB3GOs87wF+AzyrlCqPYnwizjqCI6i6JK3gMzVIoYQldIeRoJqSMEGVFGRRVpTF3IKXroG4bwdNSpFO1X0J8AINWusmrfUWoB4YBb4QpdhEnM0teOkfmSXN5Uj4/mgXI5V81kjmERQQ2rCrZZovKiJNUC8H3q+17jEvBA8o/CDwymgEJuLPbBBbU55Hmst2y4xRVR88RkQq+eJndt5D58AULqeDxprkOcYlnDnN1ywbdqMi0leheYxqveXiXsEnoqc9BSr4TBvKcklzOegbmWFuQZZN4+FM1ziBgHHkSbJs0F3OPMBQd8oIKhoiTVAPAl9QSoVqKYNffx74fTQCE/HXkQIVfKb0NCc15cY6m/nvFrFlTu8l4/qTqXFDAWkuJ10D09LrMQoiTVD/BFQAHUoprZTSQAeQA/x9tIIT8RWq4EuBBAWyYTfezAIJtTG59j+FS08zNuyCUc0n1ieiBBXsfbcbo9fdd4C7gL/WWu/TWvdFMT4RR+YxG3WVyV3BZzIr+TqkUCLmAoGlDbrmOk2y2ioHGEbNqvdBKaVeCTyktfYEvwZjzel48Gu3ed2CVkdinSamFxifXiArw0VZYZbV4cRFffDwwjZJUDE3ND7H2NQCuVnpVCf5eUlmAjZHjCJya9moex9QCQwGvz4fKZRIQOa+jdqKvKQ7pPB8QlN8fZNyEmqMha8/JfvfObzUXFpprc+qE1R4eyM7tjoS69M5YHYwT431J4CivAzyc9xMziwyND5HeVG21SElrWTeoLtcaWEWJQWZjEzM0zM0TW2SngoQD5F2knhEKfW8k8aUUmVKqUPrD0vEW1f/0ggqVTgcjlBBiKxDxVaqrD+ZzFFUc7usQ63HWtagrgO2B7+9FninUmp5P49twKbohCbiKTSCSpECCVN9dT5HW4Zp75vk8u2VVoeTlLw+Py3dEwBsSZITdC9GbSziwNFedOcYL7lio9XhJKy1rEGNAB8CHMGP9wC+sJ8HgGmMbhIiwSxN8aVWggr15JNS85jp7J9i0eOjsiSbgtwMq8OJCxlBRcda1qCOAY0ASqlHgZu01lKmkgQmZxYZn1og0+2iNEUq+EzmpmSp5IudM13Gy8SW2tSY3gNorCkgzeWgc2CK2XlP0h7+GWuR7oO6fqXkpJRyK6VeuP6wRDyFV/ClWsVRbUUeTgf0DE3j8foufgOxZme6jA2rqTK9B5CR7qKhuoBAYOnfL9YuovOglFL7gW9grEktT3KBSO9XWMNsEptKBRKmTHcaVaW59AxN0zUwTeOG5GxiaqVUTFBgrEOd6RpHd4yxZ0uZ1eEkpPUctzEEvAGYA94EfBRjDeqN0QlNxIu5/rQxxQokTEv7oSYsjiT5LHp8dPRN4nCQcslfhfZDyUpIpCJNULuBD2it7waeAwa01v+B0YfvH6IVnIiPzhQsMQ9nFkq0SaFE1LX1TuDzB6gpz0u5dZjwzuaBwEqHP4iLiTRBeQHz/+bTwCXBrx/FOAJeJJCuUIl56mzSDSeHF8ZOS4pO7wFUFGdTkOtmYnqRgdFZq8NJSJEmqKeAv1NKOYEjwCuC13cC0mM+gUzNLjI2tUCGO3V68C23URJUzJxO4QTlcDhQdWbjWJnmi0SkCeqjwFuBDwA/BLYrpVqB/wZ+EqXYRByEpvfKc1Ougs9UXpRNVkYa41MLjE8tWB1OUmnpNhLU5hRMUBDeOFb2Q0Ui0jLzg0A98P1gufle4HPA24B/jFp0IuY6U3x6D8DplJZHsTC34KXbPOK9OrUKJEzS2Xx9Ii0zvw/4J631KQidD/XVaAYm4iN8D1Qqq6/K51T7KG19k+xpkpLgaGjtmcAfgMbqfNxJesT7xWypLcThMP4WCx5f0h51HyuRTvHtBzzRDERYw9wDlWo9+JaTUvPoW+ogkZrTewDZmelsrMzH5w/Q2i3PrbWKdEPtfwI/UEr9J9CKsRcqRGt9cr2BifjoStEefMtJJV/0peoG3eXUxiLa+ybRnaNsa0je4+5jIdIE9cng55+t8DM5sDBBTM8uMjq5gDvdlfJnIW0MrsF19k/h8/lxueTIs/UyE9TmmhRPUHVF/P7JDlbHjm8AAB3sSURBVKnki0CkCaohqlEIS3QNTANQk8IVfKacrHTKi7IYHJujd3gm5dfk1mt6dpG+4RnS05yhMv5UJYUSkYu0iq9Da90B5AAvAIaD99UZvC4SQNegWWIuL8YA9VVGpZkcvbF+Znl5Y3UBaSk+Gq0pzyMnM43h8TlGJuYufgMREumJuvlKqd8Cx4GfAxUY61JHlFI1UYxvtfHsUUr9RSk1o5Q6ppS6PN4xJKKlCr5ciyOxB7NQok0KJdYtNL2X4utPYGxjMI+6b26XUdRaRPrW5vNABlDDUoHE+4Fx4ItRiGvVlFJu4NcY62GFwKeBB5VSqT2vsArdg8EpPpnOAqChWnryRcvZYMVaqq8/mczGsc2yYXdNIk1Qr8LYB9VrXtBadwLvA26IRmBrcB2QrrX+otbao7X+KXAC+Js4x5FwQiOochlBATRUm1N8MoJar1TvILHc1npZh4pEpAkql2Wl5UGuddxnpLYDp5ZdawZ2xTmOhDK/6GVwbBan00FVqSQogMqSHDLcLoYn5pmalZaSkZqaNZqjutNd8uYnSAWn+Fq6x/F4/RZHkzgiTSYPAP+qlDL75weUUmUY7Y4eikpkq5cLLG8VPAukdt30RfQMThMIQHVpDulpqb2IbXI5HaEzsaRQInJng6Onhup8KdcPys12U1uRi8frp7VHTthdrUifPe8DaoERjETwB6ATKMA4EyqeZoDlbbizMQ5PFOfRFVx/knLqc5nTfG0yzRex0AZdWX86x9aN0tl8rSItMx/QWl8JvBajOOK/gl9fFr4uFScnAbXs2tbgdXEe3cH1pxqZgjmHHF64fmaBxCZJUOeQ/VBrt+qNukqp7StcHgh+mLYppeLd6uhRwKGU+kfgLuBmjBN/74ljDAkntAdKRlDnqDdHUFJqHjEpkFjZVqnkW7O1dJI4jtHGaKWWA+Z5xg7i3OpIa72olHoF8HXgE0A7cKPWeiheMSQis4uEbNI9l9mTT1oeRUYKJM6vtiKP7Mw0hsaMDbslBal5QOharCVB2ba9kdb6OHC11XEkCp/PT9/wUpsjsSQnK53y4mwGR2fpGZpO6XOyIiEFEudnbtg9fHqI5o4xrtotCepiVp2gVmphFJz2U8CDQDnQrrUOLP89YS99IzN4fQHKirLIzIi0HWPyaqjKZ3B0lrbeSUlQayQFEhe2dWMxh08PoTvGuGp3tdXh2F5StDoSaxOa3pP1pxUtnQ0lhRJrJQUSF2Zu2G1ul3Wo1Uj4Vkdi7ZY6SEiCWomUmkdOCiQuTDbsrk0ytDoSa7RUwSfrTyuRnnyRkQKJi8vNdlNTbmzYlTdAF5cMrY7EGi3tgZIR1Eoqi3PIdLsYnZxnYnrB6nASxtnQERtSIHEh2+qNcvNTMs13UcnQ6kisgd8fCHUxlzWolTmdjtAhe7IOtXpygu7qhBJUmySoi4l2q6N84t/qSKzB8MQc84s+CnLd5Oe4rQ7HtpbWoSRBrZYUSKzO1tAIaoRAQIqeLySiGmOt9QBwpVLqeoxu4mkYHcUfkjJze+sOHfMuo6cLaQytQ8k6wWpJgcTq1JTnkpftZnRygYHRWSpLcqwOybbWtQlGa/0oRqshkSCkxdHqNGwwRlCtPZKgVkMKJFbP4XCwrb6Yp0/2c6p9VBLUBUS6DypNKfUJpdQ7w649pZT630opWR21sdD6k7yIXFB9VT5Oh1GSv+jxWR2O7UkHibXZ1iDrUKsR6TPpc8BtQGvYta8AbwE+uc6YRAyZe6DkmPcLy3SnUV2Wi88foLN/yupwbK9FjnhfE6nkW51IE9TfALdorUMVe1rrHwC3YyQpYVPdg3LMxmo1mtN8sg51UaH1p5oCiyNJDJtrC0lzOejon2RmzmN1OLYVaYI634GAYxiHFgobmpxZZGJ6kUy3i7JCaVR5MY3Vsg61WuYUn1TwrU5GuotNNYUEAnI+1IVEmqAeBf5DKVVsXlBKFQKfBh6LQlwiBsJHTw7HSqemiHCNUiixKtOzi/SPzOJOc1InU8erZk7znWwfsTgS+4o0Qb0f2AT0KKWalVKngL7gtfdFKzgRXV1SYr4mZoJq653A75fdE+dj7n9qqC6QAok1kA27FxfpPqhOpdRO4K8w9kEtAhpjH5R0QLSp0AhKevCtSkFuBiUFmYxMzNM3MsOGMvm7raQlNL0ns/trYVbyne4ck8MxzyPiMnPgn4FarfXntNb/hXGa7cekzNy+lkrMZQS1Wg2yDnVRSwUSsv60FkV5mVSV5DC/6KNNWmqtSMrMU0ioxFwq+FZtk6xDXZQ5xScdJNbOPB/qZJusQ61EysxTxILHx+DYLE6ng6pSSVCrJYUSFzY9u0jfyAzpaU7pThKBHY0lAJxslXWolUiZeYroGZwmEICqkhzS02QWdrVkL9SFLRVI5JMmayhrZiaoE63SOHYl0S4z/zfgj9EITESXbNCNTEVxNtmZaYxPLTA6OW91OLbTIvuf1mVDWS6FeRmMTy/QM7TSe/7UFs0y816gEXhvtIIT0WOWmMs0zNo4HA4plLgAM0FtkQQVEYfDcc4oSpwrogQVPN59J/B64NvA14CbgGsxjoMXNiPHvEdOCiXOTwok1m9nMEEdlwT1PBEft6G1XgTuB+5XSr0UeBtwD5AO3BWd8ES09AzKJt1IyQhqZWaBhFsKJNbFHEEdP2usQ0mXlyURJyilVD3wVoyqvRqMoolvIcnJdnz+QGh+W9ag1s4cHZjTWcJg/j0aNxRIgcQ6bKzMJycrneHxOQbH5qgozrY6JNtYU4JSSmVgTOu9DWM6z4/Re28DcI3W+ki0AxTrNzA6g8frp6Qgk+zMdKvDSTi15bm4010MjM4yNbtIXrbb6pBs4UyXnKAbDU6ngx0NJTx9sp8TrcNUFNdZHZJtrPptj1Lqq0A/8E1gCrgDqNBavxQIANIz3qbMY96lg0RkXC5n6Aj4li4ZRZlCBRKSoNYtfJpPLFnLCOpdwGmMjuW/1VrLXzJBSIn5+m2uLaS5Y4yW7nEuVeVWh2MLZrKWFkfrt3OTVPKtZC0J6nrgVuBLwP9VSh0A7sYojLCUUuoHwP8EvGGXd2utW89zk5TSETwRtq5SRlCRMl+EZR3KMDG9wODYHJluFxtkZL5umzYUkOl20Ts8w+jkPMX5mVaHZAurnuLTWv9Ra/0OoBJ4AzAKfBboCN7PTUqp/JhEeXEvAG7UWueGfUhyCursNxpR1lVa9Z8n8YUKJWSKD1haf9pUU4jLKVVn6+VyOUPHb5yQab6QNZfeaK0XtdZ3a61vwkhWfwccAO4EepVS345yjBeklMoCtgKH4/m4icLvD9AVLDGXEVTkasrzyHC7GBybY2J6wepwLCcdzKNvR3Ca71jrsMWR2EfEZeYAWutx4BvAN5RSdRgdzm+JRmDhlFJuoHiFHwUwuld4gW8ppfYDXcC/aK3vi3YciWhwbJaFRR9FeRlSfbYOLqeDxuoCTrWPcrZ7ghdsTe11qBap4Iu6nY2lABxrkQRlitrmBa11p9b601rrHdG6zzBXYpzYu/yjB8gDHscYwVVjFHH8t1JqTwziSDidwfWnjTK9t25mtdqZ7jGLI7GeOcUnFXzRozYWkZXhontwmpGJOavDsYV1jaDiRWv9GHChie4Hw76+Wyn1VuA1QMrvy+oIrT/J9N56yTqUYWRijtHJebIz06gqybE6nKSR5nKyo7GUZ04NcOTMMDdcVmt1SJZL+O3fSqlXK6VuX3bZDUjraZZGUJKg1m+pki+1Wx6F+u/VFOKUAomo2rPFmOY7cmbI4kjsISFGUBfhAr4U7Kh+COMwxSuBt1salU2EElSFTPGtV3VZLlkZLobH5xifWqAwL8PqkCxxRvY/xcyeLWWAkaCkL18SjKC01r8CPgb8P2AS+CDw18GO6ynN5w+ENunKCGr9XE4HjRtkP1Sog0SdJKho21iZT0Gum5GJeTkfiuQYQaG1/grwFavjsJuBkRkWvX5KCzLJyZIefNGwuaaQE60jtHSPc9m2CqvDibtAIMCZLqNIREZQ0ed0Oti9uYzHD/dw5Mxwyp8+kPAjKHF+Sx0kZHovWlK9UGJgdJaJ6UXyc9zSdTtGZB1qiSSoJNYpFXxRFyo1T9EEpTuM0ZPaWJTy6yOxYq5DHWsZxucPWByNtSRBJbGlPVCSoKKlqiSHnKx0RifnGR5Pvb0qunMpQYnYqCzJoaI4m+k5D609qflGyCQJKol1DsgUX7Q5nQ5UnfHibI4mUonuGAVga91KjV1EtCxV86V2VwlJUEnK5/PTHezBJ8dxR5c5emgOvlinikWPj9aeCRwOqeCLNVmHMkiCSlK9wzN4fX7Ki7LIykiKYk3bMBNUqo2gWnsm8PoC1FbkycnMMWaOoE60jjC/6L3IbycvSVBJqlMq+GLGnOJr6R7H4/VbHE38hNaf6mT9KdYKcjPYUluIx+vnaAo3j5UElaTMCj4pkIi+3Gw3NeW5eLx+2npTp+3RUgWfrD/Fw+XBfXbPnBqwOBLrSIJKUm19UmIeS6k4zRcqkJAKvrjYG5agAoHULDeXBJWEvD5/aHF1e0OJxdEkJ3MUkSoJanRynsGxObIy0qiRopu42FxTSGFuBkNjc6GK3FQjCSoJNbePMjvvpaY8l0o5DiEmzFGE7kyNSj4zETfVyRHv8eJ0OkIHYz5zMjWn+SRBJSFzznrv1tTrFRcvdRV5ZLpd9I/MMj6V/EfAm9N7sv4UX5dvD07zNUuCEkniUPMgAJdtS+1jyWPJ5XKypdZch0r+UZRU8FnjkqZynE4HJ9tGmZ7zWB1O3EmCSjLD43O0902S6Xaxo1HWn2Jpa705zZfc61A+nz/Ue7BJElRc5Wals72hGL8/wOHTg1aHE3eSoJKMOXras6WM9DSXxdEkt1RpeaQ7x1hY9FFVkpOyhzRa6bLgVP3BFFyHkgSVZA41m+tPMr0Xa03BQonTnWNJ3XX68cM9AOzfVWVxJKnJPHfs2eZB/En8PFuJJKgk4vH6OXzaKC+XAonYK8rLpLIkm/lFH209yblh1+cPcOBILwAvuqTa4mhSU11lHuXF2YxPL6Rc/0dJUEmkuX2UuQUvtRXGE1rE3q5Nyd3U80TrMGNTC1SWZMsJuhZxOBxcGRy9mm8WUoUkqCRilpen4lHkVlk6FiE5E9Tjh83R0wY5oNBCV+8xRq8Hjvam1DSfJKgk0T8yw2PPdgOy/hRPu4PHIpxoG8Xj9VkcTXR5fX7+fHQpQQnrNNUVUVaUxcjEfNIX5YSTBJUE2vsm+fBdjzM6OU9TXaGUl8dRUV4mGyvzWPT4aG5PrheOo2eGmZxZpKY8l/oq6YpvJYfDwVW7jVHUE0d7LI4mfiRBJbiTbSN85CtPMDq5wO7NpXzynVeS5pL/rPG0pyk5p/nM6j2Z3rOHq4LTfH8+kjrTfPJKloD8/gCHmgf4xHee5CNfeYKZOQ8v3FXFx9++Xw6Ss0AyrkN5vH7+crwPkOk9u1B1RZQWZjE8Mc/pJN8cbpKjVhPIyMQcjzzTxUNPd9I3PANAmsvJK6+q521/vQOXjJwssbOxBKfTwemucWbnPQn/JiEQCPCLh08zM+ehviqfWulebgvmNN+v/3SWJ470srU++fsiSoKyufkFL0+f7OfRQ9082zyAObIvK8riFS+s5yX7NsrufotlZ6bTVFtIc8cYx1tH2Le90uqQIhYIBPjefSf55WMtOBzwP/+qyeqQRJir9xgJ6sDRXt726h04k7yzvCQoG5qd9/Dc6SH+fKSXp072s7BoVIeluRzs31HJS/Zt5NKmMhkx2ciepjKaO8Y4cnooYROU1+fn6788yu+f7MDldPDBW/bK9J7NNNUVUVqQyfD4HMdahkPrn8lKEpQNBAIBugenOXJmiIOnBjh6Zhivzx/6+daNRbzo0g1ce2kNBbkyWrKjPVvK+NlDpxNuHWpkYo7Dp43n3WE9yMy8F3e6i4/efrnsp7Mhp9PBS/fX85PfN/OtXx/jix+4LqmLoiRBWcDj9dPeN4HuGKO5fYxjZ4cZnZwP/dzhgG31xVyxo5KrL9lAhXSFsL2tG4twp7vo6J9ibGqeorxMq0N6Hp8/QM/gFGe6xjnZNsrxs8P0BtcyTXWVefzdzXtkq4KN3XT9Zh4+2ElH/xT3H2jjtddsinsMgUCA2Xkvo5PzLHh8NFQXxOQgy4RLUEqpfwSu1VrfGHatDvgOsB8YBN6ntf6tRSECRqXdxMwCQ2NzDIzO0j8yQ2f/FB39k3QNTJ8zQgIozMtgz+YyLmkq4/LtFTJSSjDpaS52NpbwrB7kkYNd3HzDFkviCAQCTM4sMjIxz/D4HD1D03QNTNE9OE1b7wTzi+duJs7KcLG9oYTLtlVw2bYKOYE5AWSku3jH63bxye88xY8faObqPdWUFGRF/XH8/gDj0wsMjs0yNDZH7/A0nf1TdPZP0Ts8w6Jn6bn0jht38eoXNUY9hoRJUEqpXODjwAeBe5f9+KfAX4BXAVcDv1JKXaK1bo1lTIFAgG/+6hgdfVN4vD48Pj9z816mZheZnvMQuMBWhQ1luaiNRWzdWMT2hhLqKvNkr0mCe801jTyrB/nZHzTXX1ZLcX5ko6hAIMDo5Dz9I7OMTc0zPrXA+NQCM/MeFhZ9zC/6WPT48Hj9eLx+FjxeZua8zM57mJr1PO/NT7jyoiw21xbSVFvErs2lbNpQIGuZCWjf9kqu2FHJUyf6+e5vTvKhN+1d8fd8Pj/DE/P0j8wwMjHP5MwikzMLTM8Zz6UFz/Lnko/ZOeN5NDO3yIW2W2W6XRTlZ1JWmMXOTbEZcSdMggLuB4aAbwChvv9KqSbgMuAlWutF4BGl1L3AHcDHYhnQzLyXB/7Sjte38n/FvOx0youzKS8yPuoq89hYmUdtRV7ClyKL59u7tWLpReO+E3zwlpVfNMLNLXg52z3O2Z4JznaP0943Se/wTKgwJhI5WemUFmRSUpBFdVkONeV51JTlsrEqXyo+k8jbX7uT5/Qgf3yum/lFL3WVeVSV5DAyOU9b7wTtvZP0j86ua1Nvfo6b8qIsyoqyqSjOpq4ij7rKPGrK88jJiv1rmG0SlFLKDaxU2B/QWg8Ab9Ra9yql/pWwBAVsBzq11uGT6c3AvpgFG5Sblc7XPvxiBsdmSXe5SEtzkOlOIy/bTV52urwzTUFvf+1OntWDPHaom5fvr3/eWo7H6+dk6whHWoY41jLMma7xFc+Syst2U1WaTUlBFoV5GRTlZpCTlU6GO41Mtwt3ugt3upP0NCfudBc5melkZ6aRk5VOpts2/1uLGKosyeG2V27jO/ee4KkT/Tx1on/F3yspyKSiOJvSwiwKcjMoyHGTm5VORuh5FPwIPpeyM43XsNws61/D7PRMvhJ4dIXrPiBNa32+PvO5wOyya7NAXCoLKktyZN5ehFSW5HDz9Vv46UOar//yKB976z6Gx+foH5nhOT3EoeYBZua9od93OmBTTQGbawrZVFNIQ3U+NWW55Ga7LfxXiERx47WbeYEqp73PWNvuG56hKD+Dhup8GqoL2FCWizs9cU/Wtk2C0lo/BkSyCDMDLF8hzAam1xuTEJG4+YbNPPxMJ+19k/ztv/3heT+vq8xj79YKdm8uZXtDsUz3inWpq8ynrjI5m/naJkGtw0mgTimVpbWeC17bGrwuRNxlutN4z+v38NkfHSIrI42ywixKC7Noqivkih1VVJXKiFuI1Uj4BKW11kqpI8CnlVIfxZgqfC3wQmsjE6ls79YKfvqpV1odhhAJLeETVNDNwDcx9kANA3dorY9bG5IQQoj1SLgEpbX+1xWudQGviH80QgghYkXqoIUQQtiSJCghhBC2JAlKCCGELUmCEkIIYUuSoIQQQthSwlXxxYALoL9/5T5WQgghYifstfd5PZkkQQUbz956661WxyGEEKmsCjgbfkESFBwEXgT0YTSmFUIIET8ujOR0cPkPHIELnaonhBBCWESKJIQQQtiSJCghhBC2JAlKCCGELUmCEkIIYUuSoIQQQtiSJCghhBC2JAlKCCGELUmCEkIIYUvSSWIdlFJ7gK8Du4FW4G1a6+fthhYrU0q9DfgGsBB2+T1a6+9bFFJCUErtA+7TWpcHv3cDdwGvx+iG8gWt9WcsDNHWVvj7ZQBTwGLYr/1Za/1SK+KzM6XUS4B/B7YAg8BntdbfiNVzUBJUhIL/QX4NfBG4BrgZeFAptVFrPWlpcInjBcDntdYfsTqQRKCUcgB3AJ9b9qM7AQVsAgqAB5RSPVrrH8Q5RFu7wN9vFzCqta6Mf1SJQylVC9wN3I7x2rcX+L1Sqh24jhg8B2WKL3LXAela6y9qrT1a658CJ4C/sTashLIXOGx1EAnkTuDdwKeWXb8d+LTWekxr3Y7xAvzOOMeWCM7395Pn4erUAz/RWt+jtfYHZ4seA64iRs9BSVCR2w6cWnatGePdmLgIpZQLY2r0NqVUr1KqRSn1keC7XLGyr2ut9wLPmBeUUoUYjTZPhv2ePA9X9ry/X9ALgHKl1FGl1IBS6udKqQ0WxGdrWuvHtdbvMr9XShVjNNp+jhg9ByVBRS4XmF12bRbItiCWRFSG8ULxfaABY+763cEPsQKtde8Kl3ODn8Ofi/I8XMF5/n4AM8AB4MUY01RzwD3xiisRKaUKgHuBp4BDwctRfw7KGlTkZoCsZdeygWkLYkk4Wut+4NqwS4eVUl/GWMv7qjVRJaSZ4Ofw56I8D9dAa/2B8O+VUh8AhpRStVrrLovCsi2lVBPGGtRJ4FaWnntRfw7KCCpyJzHebYXbyrnDXHEeSqkdSqk7l112A/NWxJOotNZjQD/nPhflebgGSqlPKKW2hV1yBz/Lc3EZpdQ1GKOmXwGv11rPx/I5KCOoyD0KOJRS/4hRXnkzxpqKTA2szjjwQaVUN/Ad4FLg/cB7LY0qMf0Q+LhS6ijGlN+HgC9ZG1JC2Q1cppS6Jfj9l4D7tdZDFsZkO0qpTcB9wMe01l9e9uOYPAdlBBUhrfUi8AqMxDQKfAy4UZ7Uq6O17gFeg1HpM4lRvvpJrfUvLA0sMf0LcByjivQgxt/y65ZGlFjuAMaAFqAdYz/UbVYGZFPvAfKAzyilpsM+/g8xeg7KibpCCCFsSUZQQgghbEkSlBBCCFuSBCWEEMKWJEEJIYSwJUlQQgghbEkSlBBCCFuSjbpCxIhS6nsYXZ7P506MbtCPAnla67i0Jwo26j0AvFlrffoCv+cEngRu01rreMQmRDgZQQkRO3+P0eW5CuN4FoB9Ydc+B/w5+PXMCrePlfcDRy6UnAC01n7gE8imX2ER2agrRBwopXYCx4CG4Hk5VsWRCXQCN2itj6/yNmeBO7TWj8UyNiGWkyk+ISyklLqOsCk+pVQAeCPwUYzmm88AbwL+CaP9ziTwUa31D4O3zwM+j3FcSQB4BPj7Cxwt8QZgPDw5KaX+N/AOjCNQTgH/S2v9u7Db3IMxGnwsCv9kIVZNpviEsJ9/B/4B2A/UAc9iJKbLgV8C31BKmedAfRMjkb0M4/iSAMYx3Od78/kq4AHzG6XU64KP9SaMDtT3Az9XSuWH3eYB4K8ucJ9CxIQkKCHs5yta60e11ocxukdPY4xqNPAFjHN3GpRSjRgjolu01geDo6LbMI7mfvl57vsyjIaepnpgAegITj1+ArgJ8IT9zkmMDtVbo/KvE2KV5B2REPbTEvb1LNCutTYXi80zijKAjcGvtVLnHE2WjTGqum+F+64AhsO+/xFGpWGrUuoQximp39Vaz4X9zkjwc/ka/x1CrIuMoISwH8+y7/3n+b204O9eClwS9tEEfPc8t/EDDvOb4PEwezFGXH8G3gIcDRZ1mMzXCd+q/wVCRIEkKCES1ykgHcjRWrdorVuAPuCzGElqJf0YxRAAKKVuAt6ptX5Qa/33GCOvKeCVYbcpC7utEHEjU3xCJCittVZK3Qv8QCn1HmAI+DRGcUXzeW52CNgT9r0L+KxSagCjYnA/UBn82rSHpQP9hIgbGUEJkdhux0gmv8I4ybQAeInWevw8v38/RrUfAFrrnwMfxxh1nQY+BbxXa/1I2G2uAR7QWssUn4gr2agrRApRSmVjHGv+cq31s6v4fSfQgVEp+HiMwxPiHDKCEiKFaK1nMUZL71nlTV4LtEpyElaQBCVE6vlPYLdaVpu+XHD09DHgXXGJSohlZIpPCCGELckISgghhC1JghJCCGFLkqCEEELYkiQoIYQQtiQJSgghhC39f6ZVv4k/NCCgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = gradient(results.v)\n",
    "plot(a)\n",
    "decorate(xlabel='Time (s)',\n",
    "         ylabel='Acceleration (m/$s^2$)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can compute the maximum acceleration the jumper experiences:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "16.64087716292061 meter/second<sup>2</sup>"
      ],
      "text/latex": [
       "$16.64087716292061\\ \\frac{\\mathrm{meter}}{\\mathrm{second}^{2}}$"
      ],
      "text/plain": [
       "16.64087716292061 <Unit('meter / second ** 2')>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "max_acceleration = max(a) * m/s**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Relative to the acceleration of gravity, the jumper \"pulls\" about \"1.7 g's\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "1.6980486900939398 dimensionless"
      ],
      "text/latex": [
       "$1.6980486900939398\\ dimensionless$"
      ],
      "text/plain": [
       "1.6980486900939398 <Unit('dimensionless')>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "max_acceleration / system.g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Under the hood\n",
    "\n",
    "The gradient function in `modsim.py` adapts the NumPy function of the same name so it works with `Series` objects.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "def gradient(series, **options):\n",
      "    \"\"\"Computes the numerical derivative of a series.\n",
      "\n",
      "    If the elements of series have units, they are dropped.\n",
      "\n",
      "    series: Series object\n",
      "    options: any legal options to np.gradient\n",
      "\n",
      "    returns: Series, same subclass as series\n",
      "    \"\"\"\n",
      "    x = magnitudes(series.index)\n",
      "    y = magnitudes(series.values)\n",
      "    # units = get_units(series.values[0])\n",
      "\n",
      "    a = np.gradient(y, x, **options)\n",
      "    return series.__class__(a, series.index)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "source_code(gradient)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solving for length\n",
    "\n",
    "Assuming that `k` is fixed, let's find the length `L` that makes the minimum altitude of the jumper exactly 0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The metric we are interested in is the lowest point of the first oscillation.  For both efficiency and accuracy, it is better to stop the simulation when we reach this point, rather than run past it and the compute the minimum.\n",
    "\n",
    "Here's an event function that stops the simulation when velocity is 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def event_func(state, t, system):\n",
    "    \"\"\"Return velocity.\n",
    "    \"\"\"\n",
    "    y, v = state\n",
    "    return v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As usual, we should test it with the initial conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "0.0 meter/second"
      ],
      "text/latex": [
       "$0.0\\ \\frac{\\mathrm{meter}}{\\mathrm{second}}$"
      ],
      "text/plain": [
       "0.0 <Unit('meter / second')>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "event_func(system.init, 0, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can test it and confirm that it stops at the bottom of the jump."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd3xUVeL+8c/MpJNAIJCQkAChHTqC1CCIouuqu6z6tRfABth7QxEVFFSwrKjYC5bdRSy7KjYEpENAipRD7wm9h5D6++MO/CJSEkzmTpLn/XrNi8zNTPJkF/Nw7j33HE9hYSEiIiLBxut2ABERkWNRQYmISFBSQYmISFBSQYmISFAKcTvAn2WMCQc6ABlAvstxRESkZHxAIjDHWnuo6CfKfUHhlNMUt0OIiMif0g2YWvRARSioDICPP/6Y2rVru51FRERKIDMzk2uuuQb8v8uLqggFlQ9Qu3ZtkpOT3c4iIiKn5g+XaDRJQkREgpIKSkREgpIKSkREgpIKSkREglJAJ0kYYzoD/wQMsA0Ybq192xgTBowCLsW5UPaCtXZYILOJiEhwCdgIyhjjBb4C/mmtrQZcBYwyxrQBnsQprYY49zX1Mcb0DlQ2EREJPoEcQVUH4gGPMcYDFAJ5QA7QB+hrrd0F7DLGjAD6Ax8GItiBg7ms2LCL0BAf4WE+IsJ8RISFEBHmPA/xefF4PIGIIiIifgErKGvtDmPMKOAD4D2c5S3uwrk5KxFYUuTly4BWgcr2wifzmL0k87if93o9TlmFOsV1uMSio8KIrx5JQo0o4mtEEV89ioQaUVStEqZCExH5kwJWUP5TfNnA1cA4IA34HNjtf0lWkZdnAVGBynZOxxRy8/LJzsnnUE4+2Tl5/o/zOJSbT15+IVnZeWRl5wGHTvr1IsJ8vyus+OpRJMRFkVA9isSaVagSGVr2P5SISDkXyFN8lwBdrbUP+J9PNsa8g3N6DyCyyGujgP2BCtalVRJdWiUd9/N5+QX/v7BynCLLzsljz/4ctu3KYsuuLLbuzGKL/5GVncf6zH2sz9x3zK9Xp1Y0jevG0iSlOo3rxtIgqRphob6y+vFERMqlQBZUChB+1LE8nNl8mTiTJDb5jzfl96f8XBXi8xId6SW6mCOf/QdzjxTW1qPKa9O2/Ucek+ZuBMDn9VA/qapTWCmxNKlbneSEGHxenSYUkcorkAX1AzDMGNMPeAtoB9wM3ASsBwYbYxYC0cD9wMsBzFaqoiNDia5TjQZ1qv3hc7l5BazL2MvyDbtYsX43yzfsYsOWfazauIdVG/cwfobzuogwHw2TY2mcEkvz1Bq0aVyLqAidGhSRyiOQkyQWG2MuAYYAz+OMmh621n5ljPkeGAksxpn6/iYwOlDZAik0xEujlFgapcQ6V+GArOxcVm3ac6SwVmzYzdadWSxevYPFq3fw5eRVhPg8tGgQR/tmtenQPIE6taLd/UFERMqYp7Cw0O0Mf4oxpj6wZsKECRVqNfPd+w6xcuNu7LpdLFixDbtuJwVF/q9KrFmF9s0SaN8sgVYN4wgN0TUsESl/Nm7cSM+ePQFSrbVri36uImy3USHFxoQfKaBr/tqUvQdymGe3kr5kC/PsFjK2H+B/U1bzvymriQjz0aZxrSOvrxkbefJvICIS5FRQ5UTVKmH0aJdMj3bJ5BcUYtftJH3pFtKXbmHN5r3MWpzJrMXOvVypSVXp3DKRs9unUDuuisvJRUROjQqqHPJ5PTRPjaN5ahy9L2jO9t0HmbtsC3OWbGHBim2s2byXNZv38ukPlhYN4jinQwpprZM0yUJEyhUVVAVQMzaS8zrX57zO9cnNy2fRyh1MmreBaQszjky0GP3FIrq2TqJnhxRaNqiJV1PYRSTIqaAqmNAQH+2axtOuaTwDLsll2oLNTEjfwOLVO/g5fQM/p28gvnokZ7evy9ntU0isqVOAIhKcVFAVWFREKOd2qse5neqxefv+IwW1dddB/vWj5V8/OqcAe7ZPoWsbnQIUkeCigqokkmpGc+1fm3H1X5ry2+rtTJizgWkLNx85BfjGl4voflodLjqzIXVrV3U7roiICqqy8Xo9tG5Ui9aNatH/4lZMX7iZn+Y4pwB/nL2eH2evp32zBC7p0YiWDeO0KruIuEYFVYlFRYRyTsd6nNOxHpu37efLX1YxYfb6I9PXG6XEckmPRqS1SsTnC9jeliIigApK/JJqRXPr/7XhmvOa8u20NXw9bQ0rN+zmuTHpJNSI4h/dG3Jux7pEhOuvjIgEhn7byO9Uiw7nqvOacvFZjZiYvoEvJq8iY/sB3vxyEZ/+sIwL0lK58IxUqsdEuB1VRCo4FZQcU0RYCOenpfKXzvWZ9VsGn09aiV23i3//tJzPJ63k7PYpXHRmQ5LjY9yOKiIVlApKTsjn9ZDWOom01kksWbODzyeuZPaSTL6fuY4fZq0jrVUS157fVEUlIqVOBSXFdnh5pY1b9/Hl5FX8nO5MVZ/xWwbndqzLVX8xxFXTQrUiUjo0NUtKLDk+htsvO423Bp7DeZ3rAfD9zHX0GzaBD79dwv6DuS4nFJGKQAUlpyyuWiS3X3Yarz5wFmmtE8nJzWfshBXc/PSPfD5xJTm5+W5HFJFyTAUlf1pyfAyP9OnIiDu70aphTfYfzOW9rxfTf/gEfpq9jvyC8r0ppoi4QwUlpcbUq8HTt6TxxM2dSU2qyvbdB3n53/O5c+REZv2WQXnfvVlEAkuTJKRUeTweTm+aQNsm8fzy60bGfLeM9Zn7GPrebJrVr0GfC5vTokGc2zFFpBzQCErKhNfrocfpKYx+6GxuvqglVauEsXTtTh5+dSrPvD+bbbsOuh1RRIKcCkrKVGiIj17dGvLWwHO48lxDRJiPGYsyuPW5CXwxaSV5+QVuRxSRIKWCkoCIigjlmr825fWHepLWOpHsnHze/d9i7nlxMsvW7nQ7nogEIRWUBFTN2Ege6dORwTd1JqFGFGsz9vLAK1MYNXY++7Jy3I4nIkFEBSWuaN8sgVEPnMVlPRsT4vPw/cx13PLsBH5OX6/ZfiICqKDERRFhIfS+oDkv39uDFg3i2LM/hxc//ZVHX5/Ohi373I4nIi5TQYnr6tauyrBbu3L3lW2pWiWMRau2c+fIiXz47RKyc/LcjiciLlFBSVDweDz07FCX1x/qyXmd65GXX8jYCSu4/fmJpC/d4nY8EXGBCkqCStUqYdx+2Wk8d3s36idWZcvOLJ58eybPj0nXJAqRSkYFJUGpWWoNXrznTG74ewvCw3z8Mn8Td4yYyK92q9vRRCRAVFAStEJ8Xi7u0Yh/3tcDU686O/Zk8/ibM3jji4W6NiVSCaigJOgl1Yzm2dvO4Nrzm+Lzevh66hrueXEyKzbscjuaiJQhFZSUCz6flyvOMYy4szvJ8dFs3LqfB/45hX//aMnXckkiFZIKSsqVRimxvHRvD/7erQH5BYV89N0yHnp1Kpu37Xc7moiUMhWUlDvhoT76XdSKIf27EFctArtuF3e+MInxM9ZqFQqRCkQFJeXWaU3iGXX/WXRvW4dDOfm89tkCnnpnFrv2ZrsdTURKgQpKyrXoqDAeuLY9D1x7OlUiQ0lfuoXbnp/I9IWb3Y4mIn+SCkoqhO5tkxl1/1mc1rgW+7JyGPbBHF75z3wO5ea7HU1ETpEKSiqMmrGRPNmvC/0uakVoiJcfZq3jwVemkLnjgNvRROQUqKCkQvF6Pfy9WwOev6MbteOiWL1pD3e/OJnZizPdjiYiJaSCkgqpYXIsL97Tg04tanPgYC5D3p3Fh98u0T1TIuWICkoqrOjIUAb27UifC5vj9cDYCSt4/M0Z7N53yO1oIlIMKiip0LxeD5ee3ZghA9KIjQ5n4crt3PXCJJas2eF2NBE5CRWUVAqtG9XipXvPpFn9Guzcm83A16bx1S+rdGOvSBBTQUmlEVctkmdu7cpFZzYkv6CQt7/6jWfHpJOVnet2NBE5BhWUVCohPi839mrJw707EBkewrQFm7n3pV9Yl7nX7WgicpSQQH4zY0wi8DpwFpANvGmtHWSMCQNGAZcC+cAL1tphgcwmlUvXNknUT6rKsPdnsy5zH/e9/Au3X3YaPdolux1NRPwCPYL6CsgAEoDOQB9jzNXAk4ABGgId/Md7BzibVDJ1akUz4s7u9Dg9mUM5+Yz8eC5vfbWI/AJdlxIJBgErKGNMJ6ABcKe1NttauwboAUwE+gBPW2t3WWvXAiOA/oHKJpVXRHgI917Vjlv/rzUhPg///WU1Q9+dpetSIkEgkCOo04FFwBPGmE3GmFXAxcBBIBFYUuS1y4BWAcwmlZjH4+H8tFSG9E8jJspZcPahUVPZuivL7WgilVogC6oG0A3IxRlJXQLcD/Tyf77ob4MsICqA2URo2bAmI+7qTp1aVVibsZf7Xv6F5eu1rbyIWwJZUIeAvdbaJ6y1h6y1C4C3cU7vAUQWeW0UoC1SJeCSajrXpVo3qsnufYd45NWpTFugrTtE3BDIgloGRPln7B0WAuwCMnEmSRzWlN+f8hMJmOioMJ7s14W/dKpHTl4Bwz+cw39+Wq6bekUCLJDTzH8EtgEjjTH34RTSjcAtwGpgsDFmIRCNc+rv5QBmE/mdEJ+X2y9rQ3J8NO99vZgx45eyadt+br+sDaEhPrfjiVQKARtBWWuzgTNxrj9lAN8Bz1lrxwGPA78Bi4E5wDhgdKCyiRyLx+Ph4h6NGNi3I+FhPn5O38CgN2awZ78WmxUJBE95P21hjKkPrJkwYQLJybrJUsrGqo27GfLuLHbsySYxrgqDbuxESkKM27FEyr2NGzfSs2dPgFT/bUZHaKkjkWJomBzLyLu60zC5Ghk7DvDAK1NYsHyb27FEKjQVlEgxxVWLZPitZ9C5pbMJ4uC3ZvD9zLVuxxKpsFRQIiUQER7CI3068n9nNSK/oJBRYxfwwTdLNMNPpAyooERKyOv10PdvLbjj8tPwej189vMKRo1doDX8REqZCkrkFP2lUz0eu74jYSFefpi1jmc/nENObr7bsUQqjGLfB2WMScBZTy8eZ0uMTGCetVZ7Z0ul1aF5bZ7qn8aQd2YyY1EGT749k0ev70hURKjb0UTKvRMWlDEmBLgauBtoA+TgrPzgw1lbD2PMLOA14F/W2oIyTSsShFo0iGPYbWcw+M0ZLFy5nUdfn8YTN3ehWnS429FEyrXjnuIzxpwJLAR6A+8ATYAoa22StTYBCAPaAp8AtwPLjDE9yjyxSBBKTarGs7d3o3ZcFCs37uGhUVO0GrrIn3SiEdR9wBXW2kXH+qS1thBn9YffgNeMMW2Bp4BJpR1SpDxIrFmFZ2/vxuA3Z7A2Yy8PvjKFp/p1oW7tqm5HEymXjjuCstb2Ol45Hef1v1pr/146sUTKpxpVIxh22xk0T63Bjj3ZPPzqVOy6nW7HEimXSjJJIgpIBf5wYt1aO680Q4mUZ9GRoTzZrwvPfphO+tItPDZ6OgP7dqStiXc7mki5Uqxp5saYa4GtONek0o96zCmzdCLlVERYCI9e35EepyeTnZPPU+/MZOqCTW7HEilXinsf1DCciRINcLZnL/pIKptoIuVbiM/LPVe2o1f3BuTlF/LcmHTGT1/jdiyRcqO4p/iqAqOstevKMoxIReP1eripV0uqVQlnzPilvDZuIXuzcri8ZxM8Ho/b8USCWnFHUGOAvmWYQ6TC8ng8XH5OE269tA0eD3w0fhnvfa31+0ROprgjqOeBecaYa4C1wO9uyLXWnl3KuUQqnPO71CcmKpQRH83li0kryS8o4KZeLTWSEjmO4hbUGGA/8A2guw9FTtEZbeoQ4vPy7Idz+O8vqynIL6Tfxa1UUiLHUNyC6gB0stYuLMswIpVB55aJPNK3I8Pen8PX09aQX1jIgItb4/WqpESKKu41KAvElmUQkcqkY/PaPHZDR0JDvIyfvpZXP1tAgbbrEPmd4o6ghgHvG2NGAauA3KKftNZ+W9rBRCq605smMOiGTgx9dxY/zFpHfkEBd1zeFp9GUiJA8QvqU/+fI47xuUKc1c1FpITamngev6kzQ96dxYQ5GygoKOSuK9uppEQoZkFZa7WxoUgZadO4FoNv6sxTb89k4tyNFBTAPVe1xefTf3ZSuZ1su40SMcZournIKWjVsCZP3NyFyHAfk3/dyIiP55KXr+3VpHI70T/R7jHGjDfGnG+MOe72oMaYEGPMRcaYn3A2NhSRU9CiQRxP9UsjKiKEqQs28/xH6SopqdSOe4rPWnuRMeZiYDhQzxgzCVgMbAc8QC2cXXa7AOuBIdbaz8o8sUgF1rR+DYb0T+PxN6YzfWEGz344hwev60BoiE73SeVzwr/11tovrLVtgIuA5ThldBNwPdAeZ3XzXtbaNionkdLRpG51hg7oSnRkKDN/y2T4B3PIzct3O5ZIwBV3ksQktFOuSMA0Soll6IA0Br0xg9lLMnn6vdkM7NuRsFBNmJXKQ+cNRIJUw+RYnr4ljapVwpi7bCvPvD9bIympVFRQIkEsNakaz9zS9UhJDf8gndw8TZyQykEFJRLk6iVWZeiANGKiQpm9JFOz+6TSUEGJlAOpSdV4qn8aVSJDmbEogxc+mUe+SkoquOIudYQxJh5oDYTiTDM/QmvxiZS9RsmxPNWvC4PemM6U+ZvweT3cfZWWRZKKq1gFZYy5EXgNp5yOprX4RAKkSd3qPHFTFwa/NZ1J8zbi83m48/K22qpDKqTijqAeAN4CHrHW7ivDPCJyEs1Sa/D4jZ154u2ZTJizgRCfl1v/r41KSiqc4l6DSgFeVjmJBIeWDWsy6IZOhIV4+X7mOt74YiGFhdpPSiqW4hbUD0DPsgwiIiXTpnEtHr2hE6EhXr6dvpa3v/pNJSUVSnFP8S0AXjDG9MJZ8iin6CettQ+WdjARObl2Jp6BfTvy9Huz+O+U1fh8Xq7/W3M8Hp3uk/KvuCOoM4FZQCTOArEdijzal000ESmO9s0SeLh3B3xeD19MWsmY8Us1kpIKobhr8Z1V1kFE5NR1apnIA9e157kx6YydsIJQn5erzmvqdiyRP6Uk90ElALcDLXBGXkuBt6y1q8som4iUQNfWSdx3dTtGfjyXT36w+HxeLj+niduxRE5ZsU7xGWM64lx7uhhnP6htwN+AhcYYneITCRLd2yZz91Xt8HhgzPilfDFppduRRE5ZcUdQI4FPgVustUdObhtjRgHPAzoFKBIkzjo9hfz8Al7+93ze/d9iwkJ9XNg11e1YIiVW3EkS7YEXi5aT3ys4EyVEJIic07EeAy5pDcDozxfy46x1LicSKbniFlQGUP8YxxsAunlXJAhd2DWVG3u1AOCVsfOZPG+jy4lESqa4p/jGAG8aY+4GZvqPdQFe9H9ORILQRWc24lBuPh+NX8YLn84jNMRLWuskt2OJFEtxR1BP46wm8R9gI7AJ55rUWODRsokmIqXhinMMl/VsTEFBIc9/lE760i1uRxIpluLeB5UD3GyMuR8wwEFgpbX2YEm/oTEmFlgIPG6tfd8YEwaMAi4F8oEXrLXDSvp1ReT4rju/GTm5BXz1yyqeeX82g2/sTJsmtdyOJXJCxx1BGWMuMMaEFvn4AqArUBNn8dizihwvidFAnSLPn8QpvYY4Ey76GGN6l/BrisgJeDwebuzVgvO71Cc3r4Ah781i8eodbscSOaETjaC+BmoDW/0fH0+x94MyxvQBqgKLihzuA/S11u4CdhljRgD9gQ+L8zVFpHg8Hg8DLmnNodx8fk7fwJNvz2TogDSa1K3udjSRYzpuQVlrvcf6+FQZY1KBwUAa8J3/WCyQCCwp8tJlQKs/+/1E5I+8Xg93XtGW3LwCpszfxONvzuCZW7rSoE41t6OJ/EFxV5L42V8mRx+vZYyZW4z3+4CPgPuttZlFPhXt/zOryLEsIKo4uUSk5HxeD/de3Y5OLWpz4GAug96YzvrMvW7HEvmD446gjDE9gOb+p2cC/Y0xR9/z1Azn2tHJDAKstfbzo44f8P8ZWeRYFLC/GF9TRE5RiM/LQ73bM/S92cxbtpVBb0xn2G1nkFQz+uRvFgmQE12D2gHcD3j8j9twZtkdVohTJPcV4/tcCSQZYy7xP48BXgM6Apk4kyQ2+T/XlN+f8hORMhAa4mNg34489fZMFq7czqOvT+fZ284gvoZOYEhwONE1qEU4K0VgjJkIXOKfyFBi1trfrftvjJkPvOSfZr4fGGyMWYhzyu9+4OVT+T4iUjLhoT4eu6ETg9+cwdK1O3l09DSG33YGcdUiT/5mkTJ2omnmRf8ZdSFwyBgTdazHn8zwOPAbsBiYA4zDmYouIgEQGR7C4Js60ygllswdWTz6+jR27c12O5bICU/x7TPGJFprt+KcyjvWFp0eSjDN/DBr7WlFPs7GOX14W0m+hoiUniqRoTzVrwuPvj6NNZv38tgb03nmlq5Uiw53O5pUYieaxXc2sNP/8Vn+50c/Dh8XkXIuJiqMIf3TSEmIYX3mPh5/Ywb7s3LcjiWV2ImuQU0+1scA/uWJWgPLrbWanypSQVSLDmfogDQeeXUqqzfv4fE3ZzB0QBpREaFuR5NKqLj3QTUyxkw2xnT2X3Oa7X+sM8Z0LtOEIhJQNapG8PQtXUmoEcWKDbt54q2ZHDyU53YsqYSKu0LEKzj7Pq0FrgOScaaGvw68UCbJRMQ1NWMjefqWrtSMjWTp2p0MeWcW2TkqKQms4hZUN+Ae/yoQFwHfWGtXAG8Bp53wnSJSLiXUiOLpAWnUqBrOolXbeea92eTk5p/8jSKlpLgFlQ2EGmOq4KwqMd5/vDawpyyCiYj7kmpFM3RAV2Kjw/l1+TaGfziH3LwCt2NJJVHcgvoeZ7Q0DmetvP8ZY3r6j/23jLKJSBBISYhhyIA0YqJCmbNkC89/lE5+vkpKyl5xC6o/kI4zkrrQWnsAZ++mScDdZRNNRIJF/cSqPNU/jSoRIcxYlMELn84jv+BYt0aKlJ7i7qi7H7gLwBhT1RgTa60dXqbJRCSoNEqO5cl+XRj0xnR++XUToSFe7ry8LV6vx+1oUkEVe58nY8wtxpgNwC5ghzEmwxjzcNlFE5FgY+rVYPBNXQgP8zFhzgZe/3whhYUaSUnZKO59UPcDw3Gmm3cDugMvAg8aY+4qu3giEmxaNIhj0A2dCAvx8t2Mtbz11W8qKSkTxTrFh7NO3gBr7adFjk0zxqwDhqLVx0UqlTaNazHw+o4MfXc2/5uyGp/Xww1/b4HHo9N9UnqKe4qvFs5K40ebi3PTrohUMqc3TeCRPh0I8Xn4cvIq3vt6iUZSUqqKW1C/AZcd4/gVwLLSiyMi5UnHFrV5qLdTUl9MWsn7KikpRcU9xfc48I0xpgsww3+sC/BX4JLjvktEKrzOLRN5qHcHhn8wh88nrcTjgT4XNtfpPvnTijWCstb+APQEDuGsxXcpsBfoYK39uuziiUh5cLikfF4P4yau5INvNJKSP6+4Iyistb8Av5RhFhEpx7q0SuSh3u159sN0xk1cicfjofcFzTSSklN23ILyb6vxEs5o6RDwBfCw9n8SkePp0iqJB69rz3Nj0vns5xV4PHDd+SopOTUnOsX3JPB34DmcLTUuxFl7T0TkuNJaJ/HAde3xej2MnbCCMeOX6nSfnJITFdSlwNXW2uHW2udxZvH9wxijrTVF5IS6tnZGUodL6qPvlqmkpMROVFDJ/H4K+Rz/6xPKNJGIVAhdWyfx4LVOSf3np+V8rJKSEjpRQfmAI7uTWWsLca5FhZV1KBGpGLq2SeKBa0/H6/Xw75+W8/H3KikpvmIvFisicirOaFOH+6/xl9SPy/nke+t2JCknTjbNvK8xZv9Rr7/WGLO96Iusta+VejIRqTC6nVYHCmHEJ3P5148WjweuPq+p27EkyJ2ooNYDtxx1LBO4/qhjhYAKSkROqFvbOgCM+DidT3+w5OUXaAq6nNBxC8paWz+AOUSkEjhSUp/MZeyEFWTn5HNTr5ba9FCOSdegRCSgurWt418F3cv/pqxm1Nj52j5ejkkFJSIB17llIoNu7ERYqI8fZ6/nhU/mkpdf4HYsCTIqKBFxRTsTz1P9uhAZHsIvv25i+AdzyM3LP/kbpdJQQYmIa1o0iGPogDSiI0OZtTiTIe/MIjsnz+1YEiRUUCLiqiZ1q/PMrV2JjQ7n1+XbeOKtmWRl57odS4KACkpEXJeaVI1nbu1KXLUIFq/ewWOjp7MvK8ftWOIyFZSIBIWUhBiG33YGCTWiWLFhNwNfm8aufdluxxIXqaBEJGjUjqvC8NvOoE6taNZm7OWRV6exffdBt2OJS1RQIhJUasZGMuy2rtRPrMqmbft56NWpZO444HYscYEKSkSCTvWYCJ65tSuNU2LZujOLh1+dyoYt+9yOJQGmghKRoBQTFcbQAWm0aBDHjj3ZPPLaVFZs2OV2LAkgFZSIBK2oiFCeuLkzbZvUYs/+HAa+No30pVvcjiUBooISkaAWERbCoBs70+P0ZLJz8hny7ix+mr3O7VgSACooEQl6oSFe7r2qHZee3ZiCgkJe/vd8Pv3BanfeCk4FJSLlgsfjoc+FzRlwSWs8Hvjk+2W8+tkC8rXIbIWlghKRcuXCrqk80qcjYSFevp+5jqHvzSb7kNbvq4hUUCJS7nRplcjQAV2JiQolfekWBr4+jd37DrkdS0qZCkpEyqVmqTV47o5uxPuXRnrwlSls3r7f7VhSilRQIlJuJcfHMOKObjSoU42MHQd48JUpLF+ve6UqioAWlDHmXGPMXGPMXmPMSmNMf//xMGPMm8aYncaYbcaYRwKZS0TKr+pVIxh2a1famXjnXqnXpzF7SabbsaQUBKygjDEpwDhgKBALXAUMM8acBzwJGKAh0AHoY4zpHahsIlK+RUWEMujGTpzdPoVDOfk8/e4svp+51u1Y8icFcgRVH/jEWvuFtbbAWjsHmAR0BfoAT1trd1lr1wIjgP4BzCYi5VyIz8vdV7blinOaUFAIo8Yu4KPvllJQoHulyquQQH0ja+0UYMrh58aYGkA3YAyQCCwp8vJlQKtAZRORisHj8XDt+c2Ii41k9LgF/PvH5azP3MfdV7YlKiLU7XhSQq5MkjDGVAP+C8wC5voPZxV5SdUDqa8AAA4rSURBVBYQFehcIlIxnN+lPo/d0ImoiBBmLMrgAc3wK5cCXlDGmCbATGALcClweA39yCIviwL0t0lETlmH5rUZeVd3kuOjWZ+5j3tf+oW5y7TQbHkS6Fl83XFGTV8Cl1prs621u4BMnEkShzXl96f8RERKLDk+hpF3dadTi9ocOJjLk2/PZOyE5VrDr5wI5Cy+hsDXwOPW2kestUX/howBBhtjahpj6gP3+4+JiPwpURGhDOzbkav/YigshA+/XcpzY9K1PFI5ELBJEsBtQAzO1PJhRY6/CjwOjAQW45Tmm8DoAGYTkQrM6/Vw1XlNaVCnGiM/mcfUBZvZuHU/j17fkdpxVdyOJ8fhKe9DXf+Ia82ECRNITk52O46IBLkNW/bx9Huz2LTtANGRoTx4XXvamni3Y1VaGzdupGfPngCp/tuMjtBSRyJSqaQkxDDyrjPp0DyB/QdzeeKtGXw+caWuSwUhFZSIVDpVIkN57PpOXHGuc1Pve18vZsTHc8nO0XWpYKKCEpFKyev1cO1fmzGwbwciw3388usmHnplKlt2Zp38zRIQKigRqdS6tEpixJ3dSapZhdWb93DPi5OYtnCz27EEFZSICHVrV2Xk3c51qX1ZuQz/YA4v/WseWdm5bker1FRQIiJAdGQog27oxIBLWhMW4mXCnA3cOXISS9bscDtapaWCEhHx83g8XNg1lZfu7UHD5Gps2ZnFI69OZcz4peTlF7gdr9JRQYmIHCUlIYbn7+jOZT0bUwj856flPPDKFDZu3XfS90rpUUGJiBxDaIiX3hc0Z9itZxBfPZKVG3Zz1wuT+Xb6Gt0zFSAqKBGRE2jRII5/3ncWZ7dPISc3n9fHLeSpd2axa1+229EqPBWUiMhJVIkM5Z6r2vFQ7/ZER4aSvnQLd4yYyKzfMtyOVqGpoEREiumMNnUY9cBZnNa4Fnv25zD0vdmMGjufg1oZvUyooERESiCuWiRP9uvCzf9oSWiIl+9nruPOkROZvThT16ZKmQpKRKSEvF4Pvbo35MW7zyQ1qSqZO7IY8u4snnh7Jpu2aTPw0qKCEhE5RfUSq/Li3Wdy80UtqRIRwrxlW7n9+Z95/+vFWoWiFKigRET+BJ/PS69uDRn98Dn8pVM98gsKGTdxJbc8O4GJczfotN+foIISESkFsTHh3HH5aYy4szumbnV27j3EC5/M46FRU1m5cbfb8colFZSISClqUrc6z93RjbuvbEtsTDhL1+7k3pcm8+pnC9iz/5Db8coVFZSISCnzej307FCX0Q/15KIzG+L1ePhuxloGDJ/AN1NXk691/YpFBSUiUkaqRIZyY6+WvHL/WZzWpBb7D+Yy+otF3P3iZBas2KbrUyehghIRKWMpCTE81a8LA/t2JL5GFGsz9vLY6Ok88M8pzFi0mYICFdWxhLgdQESkMvB4PHRplUi7pvF8OXklX01ejV2/i2fen0OdWlW4uEdjzm6fTGiIz+2oQUMjKBGRAAoP9XHFOYZ3HzuXfhe1Ir56JJu2HWDU2Pnc9PSPjPt5BQcO6h4q0AhKRMQVEeEh/L1bA85Pq8/U+ZsYN3ElazP28v43S/jPhOWc36U+vbo3pEbVCLejukYFJSLiohCflx6np3Bmu2Tm2a2M+3kli1ZtZ9zElXz1y2rObp/CJWc1ok6taLejBpwKSkQkCHg8Hk5vmsDpTROw63YybuJKZv6WwQ+z1vHj7HV0bplIr24NaJ4ah9frcTtuQKigRESCjKlXg4F9O7Jx6z6+mLSKn9M3MGNRBjMWZVCjajhprZLo2iaJZqlx+CpwWamgRESCVHJ8DHdcfhpXn2f4ZtoaJs/byNZdB/l62hq+nraG6jHhpLV2yqp5BSwrFZSISJCLqxZJ7wuac935zVixYTfTFmxm6sLNbN2ZxTfT1vDNtDXExoST1iqRM9rUoXmDilFWKigRkXLC4/HQpG51mtStTt+/NWflRn9ZLdjMlp1ZfDt9Ld9OX0tsdDhdWiXStU0SLRvE4fOVzzuKVFAiIuWQx+OhcUp1GqdUp8+FzVm1aQ/TFmxm2oLNZOw4wPgZaxk/Yy1hoT5SE6vSoE41GiZXo0GdatSrXZWw0OC/IVgFJSJSznk8Hholx9IoOZbeFzRj9aY9TFu4mekLM9i0bT92/S7s+l1HXu/zekhJiDlSWg3rxJKaVJWoiFAXf4o/UkGJiFQgHo+HhsmxNEyOpfcFzdmflcPqzXtYtXEPqzftYdWm3Wzaup+1GXtZm7GXn9M3HHlvUs0qNKhTjZqxkUSFhxAZEUJk+B8fURGhzscRIYSX4UhMBSUiUoFFR4XRulEtWjeqdeRY9qE81mbsZdWmPazauJvVm/ewLmMfm7cfYPP2AyX6+ud3qc+tl7Yp7diACkpEpNKJCA+haf0aNK1f48ix3LwCNmzZx+pNe9h74BBZh/I4mJ3HwUN5zsdFnh9+ZOfkExJSdhMwVFAiIkJoiJcGdZxJFMGifM49FBGRCk8FJSIiQUkFJSIiQUkFJSIiQUkFJSIiQUkFJSIiQUkFJSIiQaki3AflA8jMzHQ7h4iIlFCR391/WDOpIhRUIsA111zjdg4RETl1icCqogcqQkHNAboBGUC+y1lERKRkfDjlNOfoT3gKCwsDH0dEROQkNElCRESCkgpKRESCkgpKRESCkgpKRESCkgpKRESCkgpKRESCkgpKRESCkgpKRESCUkVYSeJPMca0AUYDrYHVwA3W2j/c0VwZGGM6Al9ba+PdzhJIxphzgeFAY2Ar8Ly19g13UwWWMeZvwDNAKs7/Bs9Vtv8NDjPGxAILgcette+7HCegjDE3AG8Ah4ocvs1a+4EbeSr1CMoYEwZ8BfwbiAWeBn4wxlR1NViAGWM8xpibgB+AMLfzBJIxJgUYBwzF+TtwFTDMGHOeq8ECyBiTCHwGPGStjQEuA14yxrRzN5lrRgN13A7hknbASGttdJGHK+UElbyggB5AqLX2JWttrrX2X8Bi4Ap3YwXck8AtOL+kK5v6wCfW2i+stQX+0fMkoKurqQLIWpsB1LLWjjfGeIE4IA/Y526ywDPG9AGqAovczuKS04H5boc4rLIXVHNg6VHHlgGtXMjiptHW2tOBdLeDBJq1doq1dsDh58aYGjiLD//qXqrAs9buM8ZE4Zza+QF41Vq7wuVYAWWMSQUGAze4ncUNxhgfzqWO64wxm40xK40xDxtjPG5lquwFFQ1kHXUsC4hyIYtrrLWb3c4QDIwx1YD/ArNwTv1WNtlAFaADcIMx5kaX8wSM/5fzR8D91trKurlcLZx/pH6Acy3yUpwzK7e4FaiyT5I4AEQedSwK2O9CFnGRMaYJTiktAa6x1ha4HCng/D9zDpBujHkT+AfwjrupAmYQYK21n7sdxC3+Yj6zyKH5xphXgP8DXnMjU2UfQS0BzFHHmvqPSyVhjOmOM2r6ErjUWpvtcqSAMsacaYyZe9ThcGC3G3lcciVwqTFmtzFmN85p/teMMa78YnaDMaaFMebJow6H4YysXVHZR1ATAY8x5h5gFM6/FFoDX7iaSgLGGNMQ+Bp41Fr7itt5XDIfqGOMuRd4GegE3Ahc7GqqALLWNi363BgzH3ipkk0z3w3cZ4zZiDNybgvcCdzuVqBKPYKy1uYA5+MU007gUeAia+02V4NJIN0GxOBMLd9f5PGs28ECxVq7B7gAuATnv4M3gZustZNdDSYBZa3dBPQC+gN7cW6/GGKt/cytTNpRV0REglKlHkGJiEjwUkGJiEhQUkGJiEhQUkGJiEhQUkGJiEhQUkGJiEhQquw36oqcMmPM+0CfE7zkSZyV0ScCMdbagCyh5V9XbhrQ21q7/ASv8wIzgeustTYQ2URKQiMokVN3F5Dof/TwH+tY5NgIYLr/4wMBzHUnsOBE5QRH1t57Cmf/I5Ggoxt1RUqBMaYlzh5CqdbatS7miADWA2dba38r5ntWATdaayeVZTaRktIpPpEyZIzpQZFTfMaYQpxdex/BWag4HbgWeAC4DmeJmUestWP8748BRuJsfVAI/AzcdYItUq4EdhctJ2PMIKAfznYKS4GB1trxRd7zBc5ocFIp/MgipUan+EQCbzhwN9AZqAvMwymmDsDnwBvGmGj/a9/EKbLzcLZCKAS+N8Yc7x+XFwLfHX5ijLnY/72uxVmp/xtgrDGmapH3fAecc4KvKeIKFZRI4L1qrZ1orZ2Ps5L6fpxRjQVewNmjLNUY0wBnRHS1tXaOf1R0Hc429X89ztduDywu8rw+zi656/ynHp/CWRQ2t8hrluBs3vm7Fb1F3KZ/MYkE3soiH2cBa621hy8GH957Jxyo5//YGvO7bcuicEZVXx/jaycA24s8/whnpuFq/55P/wXes9YeLPKaHf4/40v4c4iUKY2gRAIv96jnx9u9N8T/2rbAaUUeTYD3jvOeAsBz+Il/65jTcUZc04G+wEL/pI7DDv8eyC/2TyASACookeC1FAgFqlhrV1prVwIZwPM4JXUsmTiTIQAwxlwC9LfW/mCtvQtn5LUPZ/+nw2oVea9I0NApPpEgZa21xpj/Ah8aY24DtgFP40yuWHact80F2hR57gOeN8ZswZkx2Bmo7f/4sDbALn5/6lHEdRpBiQS3Pjhl8iUwB6gGnGut3X2c13+DM9sPAGvtWGAwzqhrOTAUuN1a+3OR93QHvrPW6hSfBBXdqCtSgRhjooC1wF+ttfOK8XovsA5npuCUMo4nUiIaQYlUINbaLJzR0m3FfMs/gNUqJwlGKiiRiudFoLU5am760fyjp0eBAQFJJVJCOsUnIiJBSSMoEREJSiooEREJSiooEREJSiooEREJSiooEREJSv8Pk0BOE6/GAn8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "results, details = run_ode_solver(system, slope_func, events=event_func)\n",
    "plot_position(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "5.448756012469637 meter"
      ],
      "text/latex": [
       "$5.448756012469637\\ \\mathrm{meter}$"
      ],
      "text/plain": [
       "5.448756012469637 <Unit('meter')>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "min(results.y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Write an error function that takes `L` and `params` as arguments, simulates a bungee jump, and returns the lowest point.\n",
    "\n",
    "Test the error function with a guess of 25 m and confirm that the return value is about 5 meters.\n",
    "\n",
    "Use `root_scalar` with your error function to find the value of `L` that yields a perfect bungee dunk.  Hint: before calling `root_scalar`, make a version of `params` with no dimensions.\n",
    "\n",
    "Run a simulation with the result from `root_scalar` and confirm that it works."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "def error_func(L, params):\n",
    "    \"\"\"Minimum height as a function of length.\n",
    "    \n",
    "    L: length in m\n",
    "    params: Params object\n",
    "    \n",
    "    returns: height in m\n",
    "    \"\"\"\n",
    "    params = Params(params, L=L)\n",
    "    system = make_system(params)\n",
    "\n",
    "    results, details = run_ode_solver(system, slope_func, events=event_func)\n",
    "    min_height = get_last_value(results.y)\n",
    "    return min_height"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "5.448756012469637 meter"
      ],
      "text/latex": [
       "$5.448756012469637\\ \\mathrm{meter}$"
      ],
      "text/plain": [
       "5.448756012469637 <Unit('meter')>"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "guess1 = 25 * m\n",
    "error_func(guess1, params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "-1.398018644734631 meter"
      ],
      "text/latex": [
       "$-1.398018644734631\\ \\mathrm{meter}$"
      ],
      "text/plain": [
       "-1.398018644734631 <Unit('meter')>"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "guess2 = 30 * m\n",
    "error_func(guess2, params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>converged</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>root</th>\n",
       "      <td>28.996669054031372 meter</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "converged                        True\n",
       "root         28.996669054031372 meter\n",
       "dtype: object"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "res = root_bisect(error_func, [guess1, guess2], params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>success</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>message</th>\n",
       "      <td>A termination event occurred.</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "success                             True\n",
       "message    A termination event occurred.\n",
       "dtype: object"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "L = res.root\n",
    "params_solution = Params(params, L=L)\n",
    "system_solution = make_system(params_solution)\n",
    "\n",
    "results, details = run_ode_solver(system_solution, slope_func, events=event_func)\n",
    "details"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "4.1424875280200724×10<sup>-7</sup> meter<sup>2</sup>"
      ],
      "text/latex": [
       "$4.1424875280200724\\times 10^{-7}\\ \\mathrm{meter}^{2}$"
      ],
      "text/plain": [
       "4.1424875280200724e-07 <Unit('meter ** 2')>"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "min_height = get_last_value(results.y) * m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
